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SPECTRA  AND  COVARIANCES  FOR  "CLASSICAL"  NONLINEAR  SIGNAL 

PROCESSING  PROBLEMS  INVOLVING  CLASS  A  NON-GAUSSIAN  NOISE 

PART  I.  ANALYTIC  RESULTS  AND  NUMERICAL  EXAMPLES 

1.  INTRODUCTION 

Non-Gaussian  noise  fields  play  a  critical  r61e  in  modern 
signal  processing  because  of  the  frequently  dominant  effects  of 
such  noise  and  interference  in  a  wide  variety  of  applications. 
Communication  theory  generally,  and  specifically  telecommunica¬ 
tions,  electromagnetic  and  acoustic  scattering,  man-made  and 
natural  ambient  noise,  optics,  and  underwater  acoustics,  are 
common  areas  of  interest  in  this  respect.  In  the  present  report 
we  are  concerned  primarily  with  underwater  acoustic  noise 
phenomena,  but  the  models  and  results  are  canonical ,  that  is, 
they  take  forms  invariant  to  the  particular  physical  application 
in  question. 

Specifically,  we  are  concerned  with  various  second-order 
statistics  of  non-Gaussian  noise  processes  and  fields  after  they 
have  been  subjected  to  different  types  of  nonlinear  operations, 
such  as  rectification  and  modulation.  A  generic  problem  here  is 
the  passage  of  non-Gaussian  noise  through  a  zero-memory  nonlinear 
(ZMNL)  device.  The  desired  output  statistics  are  typically  the 
mean  (dc),  mean  intensity  (power),  the  covariance  or  correlation 
function,  and  the  associated  spectra.  These  last  include 
wavenumber  spectra  in  the  case  of  noise  fields,  as  well  as  the 
more  general  frequency-wavenumber  spectra  obtained  by  joint 
temporal  and  spatial  Fourier  transformations.  Typical  "class¬ 
ical"  problems  include:  (i)  rectification,  (ii)  determination  of 
output  spectra  and  covariances,  (iii)  calculation  of  ( output ) 
siqnal-to-noise  ratios,  (iv)  modulation,  (v)  demodulation,  and 
(vi)  special  systems,  as  for  example,  the  spectrum  analyzer. 

These  and  other  problems  involving  ZMNL  devices  are  described  in 
detail  in  [1;  chapters  5  and  12  -  17].  What  is  new  here  is  the 
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use  of  the  approximate  second-order  probability  density  functions 
and  characteristic  functions  in  the  above  applications  when  the 
noise  processes  are  non-Gauss ian. 

A  full  treatment  is  given  in  a  current  study  by  Middleton, 
[2],  which  is  an  expanded  version  of  his  recent  paper  [3],  which 
employs  some  of  the  results  of  the  present  report,  namely,  the 
calculated  covariances  and  spectra.  Here,  we  are  content  to 
summarize  the  pertinent  analytic  results,  the  corresponding 
examples  of  calculated  covariances  and  spectra,  and  the  various 
computational  procedures  associated  with  their  evaluation.  The 
details  of  the  derivations  are  provided  in  [2]  and  [3].  Included 
here,  also,  is  a  selection  of  illustrations  of  the  analytic 
results. 

2.  ANALYTICAL  RESULTS:  A  SUMMARY 

In  the  present  study,  we  address  three  classical  problems 
where  the  goals  are  the  calculation  of  the  covariance  and 
associated  intensity  spectrum.  Specifically,  we  consider: 

Problem  I.  The  half-wave  v-th  law  rectification  of  Class  A 
noise  fields  and  processes; 

Problem  II.  Phase  modulation  of  a  carrier  by  a  Class  A  noise 
process;  and 

Problem  III.  Frequency  modulation  of  a  carrier  by  a  Class  A 
noise  process. 

Class  A  noise,  as  noted  in  section  3  of  [2],  [3],  is  a 
canonical  form  of  interference  characterized  by  a  coherent 
structure  vis-A-vis  the  (linear)  front-end  stages  of  a  typical 
receiver:  negligible  transients  are  produced  at  the  output  of 
these  stages.  Class  B  noise,  on  the  other  hand,  is  incoherent 
and  highly  impulsive,  such  that  the  front-end  stages  of  the 
receiver  generate  an  output  which  consists  solely  of  (over¬ 
lapping)  transients.  Here,  the  Class  A  models  are  tractable  in 
the  required  second-order  distribution  and  characteristic 
functions,  whereas  the  Class  B  models  are  not  and  must 
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consequently  be  appropriately  approximated  in  second-order;  see 
[4]  and  [5]  for  additional  information.  In  the  present  report, 
we  shall  consider  examples  of  Class  A  noise  only. 

2.1  THE  SECOND-ORDER  CLASS  A  CHARACTERISTIC  FUNCTION 

In  applications  [1]  -  [3],  the  second-order  characteristic 
function,  F2 ( i^j , iJ;2 )  /  plays  a  key  r61e:  from  it,  we  may  obtain 
the  aforementioned  statistics  of  the  outputs  of  ZMNL  devices, 
spectra  of  angle-modulated  carriers,  and  other  usually  second- 
order  statistics  of  various  nonlinear  operations  arising  in  a 
variety  of  communication  and  measurement  operations. 

(See  [2],  [3]  for  further  discussion.) 

Here,  we  specifically  use  the  approximate  Class  A  noise 
characteristic  function,  F2,  including  an  additive  Gaussian 
component,  given  by 


where  A  (“A^)  is  the  "overlap”  index,  and  where 


+  T*  ^ 

A 


i]  ®2A' 


“2A  *  <2. 2b) 


^L+G  “  (”  ^G  ^a)  ®2A' 


(2.2c) 


and  k- ,  k„  are  the  normalized  covariances  of  the  non-Gauss  and 

ii  I9 

Gauss  components,  respectively.  Thus,  i  1* 
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Here,  p  (*Pj^)  is  the  "overlap"  correlation  function 


P(T') 


I 


1  -  pIt'I 

0 


for  P |t ' I  i  1 
for  P |t '  )  >  1 


(2.3) 


in  which  T  IS  the  mean  duration  of  a  typical  noise-signal  of 
.  ®  2  2 

intensity  <B  >/2  *  <L  >.  The  time  delay  t'  is  given  by 


T'  -  T  -  ^  or  T'  -  T  (»  tj-t^)  ,  (2.3a) 

o 

respectively,  for  space-time  fields  and  received  temporal 

processes.  The  path  delay  AR/c^  (*  IR2  ~  Rj^l/c^)  accounts  for 

the  time  differential  between  propagation  paths  to  the  points  at 

which  processing  occurs,  cf.  figure  2.1  ff.  Case  A.  The 
•  •  2 

quantities  22^^  and  are,  respectively,  the  intensity  of  the 
non-Gaussian  and  Gaussian  components  which  constitute  the  general 
Class  A  model  used  here.  (However,  we  note  that  the  present 
Class  A  model  belongs  either  to  the  strictly  canonical  Class  A 
cases,  where  all  interfering  sources  are  equidistant  from  the 
observer,  or  more  generally,  to  the  much  broader  class  of 
situations  in  which  the  effective  source  distribution  is 
concentrated  in  an  annulus  whose  inner-to-outer  radii  have  a 
ratio  0(1/2)  or  less.  The  former  is  exactly  represented  by  (2.1) 
to  second-order,  while  the  latter  is  approximately  so 
represented,  albeit  a  good  approximation  as  long  as  the 
aforementioned  source  annulus  is  not  too  large.  See 
[5;  section  V,  C],  for  example.  For  an  exact  treatment,  see  also 
[6],  in  an  important  class  of  physical  models.  Finally, 
differentiation  of  F2#  (2.1),  in  the  usual  way,  gives  us  the 
(exact)  covariance  of  the  composite  Class  A  and  Gauss  field, 
namely. 


K 


A+G 


‘2 


^1-52-0 


Ka  +  Kg  , 


(2.4) 


which  in  normalized  form  is 
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1  +  r' 


(2.4a) 


In  practice,  A  is  usually  less  than  unity,  say  0(0.1  -  0.3) 

mi+m2+n 

typically,  so  that  only  a  comparatively  few  terms  in  A  are 

needed  for  numerical  evaluation  of  (2.1)  and  the  statistical 
quantities  derived  from  it,  cf.  section  2.2  ff.  Note  that  when 
p I T I  1  1 ,  p  *  0 ,  and  r ' *  0 ,  we  get 


F 


2 -A 


-A 


f  S  w.3|  k  t « 


m*0 


n“0 


as  expected:  there  is  now  no  correlation  between  process  samples. 
With  a  Gaussian  component,  these  will  be  correlated,  of  course, 
unless  |t|  so  that  kg  -♦  0,  cf.  (2.2c). 

2.2  PROBLEM  I:  HALF-WAVE  v-TH  LAW  RECTIFICATION 

(STATIONARY  AND  HOMOGENEOUS  FIELDS) 


Here  we  consider  the  problem  of  obtaining  the  second-order 
(second-moment)  statistics,  M^,  of  a  sampled  noise  field,  a(R,t), 
after  passage  through  a  ZMNL  device,  g,  when  the  noise  is 
generally  non-Gaussian.  Various  processing  configurations  are 
possible.  We  show  two  in  figure  2.1,  below.  Analytically,  we 
have,  for  stationary  and  homogeneous  inputs  [1;  section  2.3-2] 


My(AR,T) 


g(Xl)g(X2) 


(2n) 


J  J  £(Uj)  £(i£2) 


=1=2 


X  F2(i^j,i52'^**'^^x  *^^2  “^1^2  ' 


(2.5) 
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where  AR  -  Rj-Rj,  r  ■  is  the  Fourier  transform 

of  the  ZMNL  device  with  ■  y(Rj/tj),  etc.  In  the  present 
cases,  we  have  specifically 

f(i^)  -  p  r(v+l)/(i^)'''^^  ,  V  >  -1  ,  (2.6) 

for  these  half-wave  v-th  law  rectifiers  [1;  (2. 101a, b)]. 


SPATIAL  SAHPLING  -*•  ^  SIGNAL  PROCESSING  - ► 


2 -SENSOR  ARRAY  ZMNL 


Figure  2.1  A.  Two-point  sensor  array  (Rj)  giving  sampled  field 
at  two  space-time  points.  B.  A  general  array  (R)  (preformed 
beam),  converting  the  field  a(R,t)  into  a  single  (time)  process 
x(tj).  Both  are  followed  by  ZMNL  devices,  delays,  and  averaging, 
as  indicated  schematically. 
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For  the  Class  A  non-Gaussian  noise  inputs  of  section  2.1 
above,  we  find  that  the  (normalized)  second-moment  for  the 
resulting  rectified  field  is  now 


mi+m2  » 


My(AR,T)  -  exp[-A(2-p)] 


,  [A(l-p)] 

^  Z—  n*!'  ®2* 

mi,m2*0  n*0 


+  r' 


v|m^,m2/n  ' 


(2.7) 


where  we  have  further  postulated  the  noise  field  to  be  isotropic, 
AR  -»  |AR|,  and  where  specifically. 


(2.7a) 


fm^+n 


+  r' 


m2+n 


+  r' 


;  a  »  (mj,m2,n)  ,  lY^^I  i  1  . 

(2.7b) 


Specifically,  also,  we  have  the  following  normalized  forms 
My  s  My/2^^2V4n  ;  f '  «  pr '  ,  p  -  l/T^  ,  cf.  (2.3)  , 


^  m  AR/Aj^  ,  Aj^  ■  correlation  distance,  AR  *  |R2-R]^1  .  (2.8) 


For  numerical  results,  we  select  the  following  models  for  the 
space-time  covariance  functions  of  the  isotropic  and  stationary 
non-Gaussian  and  Gaussian  components  of  the  input  noise  field: 


*  exp(-AR^/Aj  -  j(AWj^T'/p)^)  «  exp(-^^  -  i(AUj^f')^)  ;  (2.9a) 
kg  *  expl-AR^/A^  -  |(AWgT'/p)^)  -  exp(-AR^(Aj^/Ag)^  -  ^(AUgf')^]  ; 
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,  AWg  *  AcOg/p  .  (2.9b) 

Here,  Ag  is  a  correlation  distance,  and  Aw^^,  Aw^  are  angular 
frequency  spreads  associated  with  the  respective  non-Gaussian  and 
Gaussian  components  of  the  input  field.  Note  that  if  we  define 
the  correlation  distance  A^^  as  that  where  *  1/e  (f'  »  0), 
then  Aj^  ■  etc. 

For  the  special  cases  of  v  considered  here,  we  also  observe 
(from  [1;  (A. 1-39)])  that  may  be  expressed  in  closed  form; 


BQ(y)  =  n  +  2  arcsin(Y) 

9 

(2.10a) 

Bj(Y)  *  Y  arcsin(Y)  +  |l  -  Y^ 

+  ^^Y 
“  2  ^  ^ 

(2.10b) 

B2(Y)  '  (i  +  +  arcsin(Y))  + 

|y(i  -  . 

(2.10c) 

2.2-1  GAUSS  PROCESSES  ALONE  (A»0) 


When  only  a  Gauss  noise  field  is  originally  present,  that  is, 
A  *  0,  for  example,  S. 
result  [1;  page  541,  (13.4a)]: 


0,  for  example,  ®2A  “  (2.7)  reduces  to  the  classical 


(2.11) 


For  comparison  with  the  non-Gaussian  cases  (A>0),  we  choose  to 
have  equal  input  noise  intensities.  This  means  that 

Vo  -  ”1  *  V  -  V'  r'l  ' 

so  that 


“yIa-O 


(1  +  r')'’ 


a*0 


(2.12) 


A  ^ 

and  My  is  then  to  be  compared  with  My,  A  >  0.  When  r'  is  small, 
as  is  usually  the  case,  we  can  often  replace  (1  +  r')'*  by  unity. 
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At  this  point,  following  figure  2.1,  we  distinguish  two 
classes  of  operation:  (A),  where  a  pair  of  point  sensors  is  used 
to  sample  the  noise  field  and  we  wish  t  consider  both  the  space 
and  temporal  correlations  of  the  sampled  field  at  the  two  points 
(Rl,ti),  (R2,t2)?  and  (B),  where  the  space-time  field  is 
converted  into  a  random  process,  x(t),  by  the  beamforming  array 
(R),  with  an  associated  directionality  embodied  in  the  resultant 
beam  (vide  [7;  sections  IV  B  and  VI  A]). 


2.2-2  CASE  B,  FIGURE  2.1 

Let  us  consider  the  simpler  case  (Case  B)  of  the  time  process 
first,  cf.  (B).  For  this,  we  set  AR  *  0  formally  in  (2.7)  et 

A 

seq.  above,  since  x  *  R  a(R,t)  here  and  t'  *  r  *  t2-tj^,  cf. 

(2.3a).  See  also  [3;  (3.2)  et  seq.  and  (3.11a)].  Then  our  ad 

hoc  illustrative  models  of  the  process  covariances  k^,  k^,  are, 
from  (2.9a,b),  at  once* 

\  “  exp(-  *  exp(-  -(AUj^f)^)  ,  (2.13a) 

kg  »  kg(T)  -  exp[-  ^(AWgf/p)^j  -  exp(-  |(AUgf)^]  .  (2.13b) 

Accordingly,  (2.7)  reduces  to 

Case  B:  My(0,f)  *  My(f)g  -  (2.7),  with  «  (2.7b), 

and  (2.13a,b)  and  ^  «  0  therein.  (2.14) 

We  note  that  when  |f|  21,  p  *  0,  and  My(0,|f|  2  1)  reduces  to  a 
simpler  relation  [vis-A-vis  (2.7)],  viz.; 


*  A  physically  derived  model  of  kg  and  kj^  may  be  made  from 
[3;  (3.11a)]  with  L-RL,  R-  (2.9)  etc.,  where  L  is  typically 
given  by  [3;  (3.3)],  for  example. 
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.-2A 


^  1  2  m. 

e  >  — -  — =•  +  r' 

l_ _  m.  !  m,!  A 


v/2 


iHm 

-4  +  r' 

A 


v/2 


®vla  ' 


where  (2.7b)  becomes 


(2.14a) 


r"  k. 


m. 

T~ 

m. 

— i  +  r' 

-4  T' 

A 

1  A  ) 

p  «  u 


in  B. 


V  a 


Special  cases  of  interest  are: 

I.  THE  INTENSITY  E(y^):  f  *  0,  p  -  1, 
(2.14)  reduce  to 


(2.14b) 


m2  *  0 ,  and  (2.7), 


y  norm  ■  "  “ylO'O) 


B  e 

V  a*n 


-A 


n-O 


where  now  Y^  „  ■  1,  e.g.,  k, (0) 
a*n  ^  L '  ' 

of  n,  for  example,  for  Y.  «  1, 

A 


(2.15) 

1  etc.,  and  B^  is  independent 


B 


2n 

for 

V 

«  0 

1  *  < 

n 

for 

V 

«  1 

1  a*n 

,3n/2 

for 

V 

*  2, 

,  cf.  (2.10) 


(2.16a) 


For  general  v,  Y  ■  1,  we  have  (from  [1;  (A.1-34)]) 


®v  a-n  “  2n'’r(v+J5)  ,  v  2  0  . 


(2.16b) 


Thus,  (2.15)  becomes 


y  norm  '  «yf®>B  *  «y«>'‘>l 


2«‘>r<v+i,)  (2  +  r')'’  . 

(2.17) 


n»0 
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The  unnormalized  form  is,  from  (2.8), 


My(0)B  -  My{0,0) 


^2A  (2.18) 


with 


H<''l(A,r')  IE  (2  +  r') 


1  +  r' 


for  V  »  0  , 
for  V  *  1  , 


1/A  +  (1+r')  for  V  «  2  . 


(2.18a) 

(2.18b) 

(2.18c) 


For  other  values  of  (>0),  we  must  evaluate  numerically. 


II.  THE  MEAN  VALUE,  y;  |f  |  « 

Now  p*0,  n«0,  y^*0,  and  (2.7)  reduces  directly,  upon 
use  of  (2.18),  to 


(2.19) 


The  unnormalized  form  of  (2.19)  is,  from  (2.8), 


y  -  My(»)B  -  My(0,«) 


^  “2A 


r^(^)  r  (2.20) 


and  for  v  even,  we  find,  from  (2.18a,b,c) 


-  1  ,  -  1  +  r"  ,  “  5  •  (2.21) 
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2  — 

III.  THE  CONTINUUM  INTENSITY:  y  -  y 

From  (2.18)  and  (2.20)  we  get  at  once  the  general  result 
for  V  i  0, 


which  is  the  generalization  of  [1;  (13.7)],  in  the  classical 

purely  Gaussian  cases,  to  the  present,  dominant  non-Gauss ian 

2 

noise  component  22a  these  classical  cases,  we  can 

show  at  once  that 


t.(^) 

“l 


(2.23) 


.  2 

where  22^^  •  implies  A  0  and  Bq  0,  cf.  (2.2b),  so  that 

(2.22)  becomes,  as  expected, 

P  L  -  2'’  o?'’  - 

c  I  Gauss  G  2n** 

Figure  13.5  of  [1]  shows  (2.24)  as  a  function  of  rectifier 
law  (V),  as  well  as  (2.18),  (2.20)  in  these  Gaussian  cases.  In 
the  present,  more  general,  situation  of  Class  A  noise,  the 
results  are  more  complex,  as  expected,  with  now  two  additional 
parameters  (A,r'),  descriptive  of  this  much  broader  class  of 
interference. 
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2.2-3;  CASE  A,  TIGURE  2.1 

We  turn  now  to  the  more  general  problem  of  the  covariance  of 
the  Class  A  non-Gaussian  random  field/  sampled  according  to 
procedure  (A),  shown  schematically  in  figure  2.1  earlier.  Here, 

X  ■  a(R/t),  sensed  at  (R^/tj),  (R2/t2)/  where  L  *  L,  cf.  (3.3)  in 
[3;  (3.2)].  Equation  (2.7)  applies  here,  with  AR  0  (as  well  as 
for  AR  »  0),  and  we  use  (2.9a,b)  for  our  illustrative  examples, 
which  are  discussed  in  section  3  following.  At  this  point,  we 
recall  from  (2.3a)  that  the  proper  time  delay  to  use  is 
T'  »  T  -  AR/Cq  in  p  «  p(T'),  and  in  some  of  the  structural 
elements  of  the  noise  field  covariances,  cf.  [3;  (3.11b,c)]. 


CASE  I;  f'  -  0 


From  (2.7),  we  have  p  *  1,  m^^  *  m2  “  1,  giving 


My(A'R,0) 


n-0  ° 


25) 


where  (2.7b)  is  specifically 


n 


j  kj^(AR,0)  +  r'  kg(AR,0) 


a»n 


5  +  r" 

A 


(2.25a) 


For  calculations,  (2.9a,b)  are  used,  with  given  by  (2.7a), 
where  (2.25a)  provides  Y^.  When  AR  *  0,  (2.25)  reduces  to  (2.15) 
et  seq.  for  the  total  intensity  of  the  field  observed  at  Rj^  -  R2« 
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CASE  II:  AR  -»  CO,  |t  '  I  >  1 

When  AR  -»  00,  we  obtain  different  results,  depending  on  t  '  . 
Here  p  *  0,  0,  cf.  (2.9a,b)  in  (2.7b),  and  therefore  n  «  0. 

Accordingly,  (2.7)  becomes 

2 

My(<»,|T'!  >  1)  «  My(co,«,)  .  My(0,->)  -  ,  (2.19)  .  (2.26) 

The  fact  that  AR  ->  oo  ensures  that  Y  -»  0,  a  behavior  similar  to 

a 

that  for  Case  (B)  above,  when  we  consider  the  purely  Gaussian 
noise  process,  section  2.2-1. 

CASE  III:  ^  CO,  0  <  |t  '  1  <  1 

Here,  p  >  0  while  Y  -»  0,  so  that  -  (2.7b),  becomes  r^(v+J5) 

O  V 

once  more.  The  second  moment  function  (2.7)  is  now 


(2.27) 


which  is  a  minor  simplification  of  (2.7). 

CASE  IV:  ^  CO,  |t  '  1  »  0 

In  this  special  situation,  where  x  ■  co  in  such  a  way 

that  T'  ■  0  and  therefore  p  ■  1,  Y.  -  0,  we  obtain  directly  from 

A 

(2.27)  the  comparatively  simple  result. 


My(»,0) 


(>  0)  . 


(2.28) 
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2.2-4:  REMARKS 

At  first  glance,  as  AR  -»  ®,  we  might  expect  M  always  to 

—2  —2  ^ 
reduce  to  y  ,  e.g. ,  s  -  y  »  0  for  the  covariance  of  the 

rectified  space-time  field.  This  is  expectedly  the  case  for  the 

covariance  (and  s econd -moment )  function  of  the  input  Class  A  and 

Gauss  noise  field  components  a(R,t),  as  we  can  see  directly  from 

(2.9a,b),  or  from  [3;  (3.11b,c)]  for  example,  in  the  physically 

derived  cases.  However,  the  process  or  field  y  *  g(x)  here  is 

the  result  of  a  nonlinear  operation,  cf.  (2.5),  (2.6),  which 

severely  distorts  the  input  waveform  and  generates  all  kinds  of 

modulation  products,  associated  with  the  spatial  as  well  as  the 

temporal  variations  of  the  input  field.  This  accounts  for  the 

_2 

departures  in  Cases  III,  IV  of  M  (®,f')  from  y  ,  while  certainly 

—2  —  ^  . 

Mx(®,T')^x  *0,  (since  x  =  0  initially  here). 

From  the  various  limiting  results  above,  we  see  that 

My(0,0)  >  My(«,0)  and  My(0,0)  >  My(0,®)  ,  (2.29a) 

and 

^y(®»0)  =  My(0,®)  depending  on  A,  r',  and  v  ,  (2.29b) 

with 

My(0,0) 

whereas 

_2 

M^(0,0)  >  M^(0,*)  -  M^(®,0)  -  X  «  0  .  (2.29e) 

Finally,  we  note  that  (2.11),  (2.12)  apply  here,  also,  for 
the  Gauss-alone  cases,  where  now 

2j^(i  +  r')  -  and  .0  ■»  (1  +  r')'’  -  0  . 


-  My(0,«)  >  0  ,  cf.  (2.19)  and  (2,20)  ,  (2.29c) 

2 

-  My(0,®)  -  y  ,  cf.  (2.20)  and  (2.26)  ,  (2.29d) 
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2.2-5:  SPECTRA 

The  various  intensity  spectra  associated  with  the  output  of 
the  processor  (cf.  figure  2.1)  are  important  also,  as  they  show 
how  the  energy  in  this  output  is  distributed.  Here,  we  consider 
two  types  of  spectra,  respectively,  for  the  rectified  spatial 
field  (A)  and  for  the  process  (B),  namely  the  wavenumber  and  the 
frequency  spectrum  of  y(R,0)  and  y(0,t).  In  particular, 
wavenumber  spectra  are  useful  in  the  analysis  of  spatially 
distributed  phenomena,  paralleling  the  analysis  of  time -dependent 
phenomena . 


I :  WAVENUMBER  SPECTRUM 


The  wavenumber  intensity  spectrum  is  defined  here  by 

W2(k,0)y  -  W2(k,T)y|^^Q  S  JJ  My(AR,0)  exp{ik-AR)  d(AR) 

AR 

«0 

-  2n  r  M„(AR,0)  J  (kAR)  AR  d(AR)  «  W,(k,0)  , 

j  y  o  z  y 

0 


with 


(2.30a) 


(2.30b) 


k  -  “  |k| 


(2.30c) 


for  these  isotropic  fields,  where  k  is  an  (angular)  vector 
wavenumber.  Using  the  normalization  of  (2.8),  we  get,  with 
k  s  kAj^, 


W2(k,0)y  s 


W2(k,0)y 

A?  2'*  2^'’/4n 

ii 


2ii  I  My(x,0)  Jjj()cx)  X  dx 
0 


(2.31) 


for  the  normalized  wavenumber  intensity  spectrum. 
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A 

Since  My(®,0)  IS  nonvanishing,  cf.  (2.28),  there  is  a  dc 
component,  or  6-function,  in  the  general  wavenumber  spectrum. 
We  will  use  the  relations 


PA  1  A  A  /AO  AO\^  .A 

J  X  J^{kx)  dx  -  i  5(k  -  0)  ,  k  -  (k;;  +  k‘)  -  |kj 
0  ^ 


A  A 

k  -  (k,^)  , 
(2.32) 


where  we  must  remember  that  k  is  two-dimensional.  With  v  a 
vector  wavenumber  defined  by 

A  A 

k  -  2nv  [-  (<),♦)]  ,  k  »  2nv  «  2n|v|  ,  (2.33a) 


and  using  the  relation 

6(ax  -  b)  -  ^  6(x  -  I]  for  a  >  0  ,  (2.33b) 

we  also  show  that 

*  y  2nk 

- ^  •  |v| 

(2n)'^  *  y 

Applying  (2.32)  -  (2.34),  with 

00 

W2(k,4)y  -  2n  J  x  J^{kx)  [My(x,0)  -  My(«,0)]  dx 
0 


-  - 6(v-0) 

(2n)^v 

/  V  “  (vx  Vy)**  •  (2.34) 


+  2n  My(<»,0)  J  x  JqC^x)  dx 
0 


(2.35a) 


-  W2(k/0)y_^Q^^  +  (2n)2  MyC,0)  6()^^-0)  6(}^y-0) 

-  ' 


(2.35b) 
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which  defines  the  continuous  portion  of  the 

spectrum  and  shows  the  dc  term  in  k-  or  v-space,  as  convenient. 

with  which  we  are  concerned  in  the  specific 


It  is  W 


2-cont 

numerical  examples  of  section  3  ff. 


II.  FREQUENCY  SPECTRUM 


Here  we  employ  the  Wiener-Khintchine  theorem  [1;  (3.42)]  to 
write  for  the  frequency  spectrum  of  y 

e» 

Wy(f)  “  2  J  My(0,r)  exp(-iwT)  dr  -  j  My(0,f)  cos(uf)  df,(2.36) 
—00  0 


where 


®0  *  “2A  * 


w  *  2nf;  u 


w/p;  f  *  f/p.  (2.36a) 


Accordingly,  we  define  the  normalized  frequency  intensity 
spectrum  of  y  as 


Wy(f)  *  Wy(f)/BQ  «  J  My{0,f)  cos(«f)  df  . 

0 


(2.37) 


A  _2 

Again,  there  is  a  dc  component,  since  My(0,®)  «  y  (>  0), 
cf .  (2.20) .  We  have 


Wy(f)  «  J  [My(0,f )-My(0,®)]  c o s ( )  dt  +^My(0,®)  ^(f-O),  (2.38) 
0 

since 

m 

J  cos(ux)  dx  ■  n  6(w-0)  ■  ^  6(f-0)  . 

0 
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As  in  the  wavenumber  cases  above  (Case  I)/  we  are  concerned  with 
the  continuous  part  of  the  spectrum,  viz. 


^(^^cont 


J  [My(0,f)  -  y  ]  cos(wf)  df  , 
0 


(2.39) 


which  is  also  illustrated  numerically  in  section  3  ff. 

III.  WAVENUMBER  FREQUENCY  SPECTRUM 

The  wavenumber  frequency  spectrum  is  defined  by 

CD 

W2(k,w)y  *  JJ  My(AR,T)  exp(ik*AR-iwT )  d(AR)  dr  ,  (2.40) 

— «0 

with  (0  *  2nf.  The  associated  wavenumber  spectrum  W2(k,0)  used  in 
(2.30)  is  obtained  from  W2(k,T)j^^Q.  In  normalized  form,  we  have 
for  (2.40),  in  these  isotropic  cases, 

W2(k,w)y  -  (2'’  fi2A  ^  W2(k,«)y 

00 

-  JJ  My(AR,f)  exp(ik*AR-iwf )  d(AR)  df  ; 


W2(k,w)y 


00  00 

0  0 


A  A  .  . 

My(x,f)  JgCkx)  exp(-i«T)  X  dx  dr  . 


(2.41) 


The  various  dc  components  are  readily  extracted,  as  in  Cases  1 
and  II  above.  Numerical  examples  of  this  joint  intensity 
spectrum  are  reserved  to  a  possible  subsequent  study.  The 
results  of  section  3  show  the  marginal  spectra  (Cases  I,  II)  of 
this  more  general  situation. 
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2.2-6:  FREQUENCY  AND  PHASE  MODULATION 

BY  CLASS  A  AND  GAUSSIAN  NOISE 

This  is  a  Case  (B)  situation,  cf.  figure  2.1,  where  AR  *  0 
and  we  are  concerned  only  with  the  received  (non-Gaussian)  noise 
process  which  is  used  to  angle-modulate  a  (high  frequency) 
carrier  f^.  For  the  analysis,  see  [3;  section  II]. 

The  general  result  for  the  covariance  of  the  carrier 
modulated  by  Class  A  and  Gauss  noise  is  found  to  be 


(2.42) 

where  now,  cf.  (1;  ( 4 .2 ) , ( 14 . 14c) ] , 


2(T)g  . 

’g 

i“2A) 

J  (It|-X) 

k(X)  dX 

°f) 

AjFM 

0 

-  (' 

or 

CO 

ill 

0 

'^(G  or 

L)(f)  — 

cos ( WT ) 

df  . 

(2. 

42a) 

Also, 

cf.  (1; 

(14, 

•  2  ) ,  1 

(14.14c; 

)], 

2(^)g 

or 

a“2A, 

I  {k(0) 

or  L) 

(Do  -  Dp) 

(2. 

42b) 

A 

PM 

and 

o  o 

|fm 

00 

iJ^(f)  df/w^  or 

Q  1  ae 

“o|PM 

®2A  • 

(2. 

42c) 

0 
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For  our  numerical  examples,  we  use  the  RC-spectrum  of 
[1;  section  14.1-3],  where  now 


kj^(C)  -  exp(-b|C|  ) 
and  therefore 


kg(C)  •  exp(-|C,|)  ,  C  “  T  AUjj  ,  (2.43) 


FM: 


PM: 


D|a(T)j.- 

("f) 

^  [exp(-b|C| )  +  b|C|  -  1] 

Dp  “(^>G  -  : 

-("^Ia 

[exp(-KI)  +  Id  -  1]  ; 

1''f)a  “  ®2A 

• 

i 

D?  a(T)^  - 

[1  -  exp(-b|C| )]  ; 

2(T)g  - 

*■'  ("pIa 

[1  -  exp(-b|d)]  ? 

(4)  A  -  ''P 

®2A  ' 

(2.44a) 


(2.44b) 


with  b  (>  0)  a  dimensionless  quantity,  as  is  The  quantity 
is  the  bandwidth  of  the  modulating  (Gauss)  noise,  cf. 
(2.43).  Note,  also,  that 


T' 


”1  -  (''?)g 


■  "g  ■  ("IIg  •  <2-«) 


2 

The  quantities  (Pp  p)()  respective  modulation  indexes  for 

FM  and  PM,  cf.  [1;  chapter  14]. 

Finally,  we  have  for  p  in  (2.3),  now  with  AR  ■  0, 


p(T)  ->  p(C) 


for 


A(t), 


S  1 


N 


(2.46) 
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Putting  the  above  (2.43),  (2.44)  in  (2.42),  we  now  specialize 
our  results, 

“  i^o  ^o<^)  cos(a>^T)  ,  with  k^{0)  -  1  ,  (2.47) 

to  the  normalized  covariance  ^^(t),  respectively,  for  FM  and  PM, 
and  their  associated  spectra.  We  have  for  these  carriers 
modulated  by  a  sum  of  Gaussian  and  Class  A  noise: 


I .  FREQUENCY  MODULATION 


with  p(T)  given  by  (2.46).  Here,  2^  •*  •*>  in  (2.42).  Since 


lim  k  (C) 


FM 


0  , 


there  is  no  dc  in  and  hence  all  the  original  carrier  power 

(~A^/2)  is  distributed  into  the  sideband  continuum  for  this 
'  o  ' 

highly  nonlinear  modulation,  as  expected  [1;  section  14.1-2]. 

The  associated  intensity  spectnim  for  is  defined  by 


W(u)) 


J  ’'o^^^FM  cos(wC)  dC.  , 


(2.49) 


which  is  determined  by  a  direct  cosine  transform  of  kQ(C)pjj*  See 
appendix  A. 6  ff. 
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II.  PHASE  MODULATION 


“  ®*p["^'(^p)a  +2A(1-p)  exp(-  ^[^p];^] 

-  A(2-p)  +  Ap  exp(-  [1  -  exp(-blCl)]]]  , 


as  before;  that  is,  the  total  (normalized)  intensity  is  unity. 
Also 


(2.51b) 


this  is  the  fraction  of  the  power  remaining  in  the  carrier,  so 
that 

•'o^^^PM  ■  ^o^*)pM  *  ^  "  (2.51b)  ,  (2.51c) 

which  is  the  fraction  of  the  power  distributed  in  the  sideband 
continuum. 

The  associated  intensity  spectrum  of  the  sideband  continuum 
is  determined  from 


'^^“U+GlPM-cont  "  J  [^o^^^PM  "  ^o^*^Pm]  •  (2.52) 

0 

See  seccion  3  ff.  for  examples  and  appendix  A. 5  for  the 
evaluation  methods. 
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Finally,  in  the  equivalent  Gaussian  cases  (Gauss  noise 
modulation  of  equal  intensity  and  basic  spectrum,  e.g., 

-  r'  4^  (1  +  r')  and  Kg  -  kj,  , 
we  see  that  (2.48),  (2.50)  reduce  to 

''o<'^>FM-Gauss  -  ^pI-I^^fIa  p'  texp(-bHI)  +  b|Cl  -  D/b^]  , 

(pf)a  -  (1  *  p''(''f)a  >  (2-“®' 

''o<^>PM-Gauss  ■  t*  -  eKP(-b|CI)j]  . 

(p^Ia  .  (1  +  r')(p^)A  .  (2-53b) 

with  spectra  obtained  as  before,  from  (2.49)  and  (2.52). 
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3.  NUMERICAL  ILLUSTRATIONS  AND  DISCUSSION 


It  is  convenient  to  discuss  the  general  results,  namely  the 
effects  of  (ZMNL)  nonlinear  rectifiers  on,  and  modulation  by,  a 
mixture  of  Gaussian  and  non-Gaussian  noise  processes  and  fields, 
from  the  specific  numerical  calculations  presented  here  in 
figures  3.1  -  3.10.  These  constitute  a  representative  selection 
from  the  universe  of  possible  parameter  states  [cf.  "Summary  of 
Normalized  Parameters"  and  section  2,  preceding].  This  is  done 
here  on  a  per-figure  basis,  as  noted  below.  In  each  case,  the 
dc  component  is  removed:  only  the  covariance  or  continuous 
spectrum  is  calculated.  We  recall  that  there  are  two  cases  to 
distinguish:  Case  A,  t'*t-AR/Cq,  a  2-element  array;  and  Case  Z, 
T'*T«t2-t2,  a  preformed  beam.  See  figure  2.1  and  (2.3a). 

All  spectra  shown  here  are  normalized  to  have  area  (under  the 
spectrum  level)  of  unity,  i.e.,  the  spectral  normalization  is 
obtained  by  dividing  the  spectrum  by  the  value  of  the  associated 
covariance  at  its  origin.  The  normalization  of  the  covariances 
themselves  is  obtained  by  dividing  by  the  value  at  f*0  or  AR=0. 

I.  GAUSS  NOISE  ALONE 

FIGURE  3.1 

This  figure  shows  the  normalized  temporal  covariance  (AR*0) 
for  both  the  input  and  output  of  a  ZMNL  half-wave  v-th  law  (viO) 
detector,  when  v  «  0,1,2  and  when  only  Gaussian  noise  (A«0)  is 
applied  to  these  nonlinear  devices.  These  curves  are  based  on 
(2.11)  with  (2.7a),  where  »  k^,  (2.96),  with  s  A«g/p  *  5 
here.  The  normalization  is  with  respect  to  the  covariance 
maximum;  e.g.,  the  normalized  covariance  shown  in  figure  3.1  is 
obtained  from  [ (2. 11 )/(2. 11 )f«0] ,  AR»0.  These  results  apply  for 
both  cases  A,B  of  figure  2.1,  where  now  f'»f,  since  AR»0,  cf. 
(2.3a)  and  remarks . 
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As  expected  (cf.  [1;  chapter  13]),  the  general  nonlinearity 
(2.6),  viO,  contracts  the  covariance,  which  is  equivalent  to 
spreading  the  spectrum  vis-A-vis  the  input  ,  cf.  figure  3.2, 
below.  Moreover,  the  greater  the  distortion  (v*0,2),  usually  the 
greater  are  these  effects.  [See  appendix  A.I.] 

FIGURE  3.2 

This  is  the  same  situation  as  shown  in  figure  3.1,  except 
that  the  normalized  intensity  (frequency)  spectrum  is  calculated 
now  [cf.  section  2.2-5,  Case  II,  (2.39)]  with  My(0,f),  (2.11), 
used  in  (2.39).  Observe  the  greatly  broadened  spectra, 
particularly  at  the  low  spectral  levels,  where  the  greater  spread 
occurs  for  the  "super-clipper",  v=0,  cf.  remarks,  figure  3.1; 
also,  appendix  A. 3. 

FIGURE  3.3 

For  the  same  purely  Gaussian  field  above,  cf.  (2.11)  and 
(2.96),  with  the  spatial  covariance  is  calculated,  with 

parameters  -  5^,  using  (2.11)  as  before.  The  normalization 

is  with  respect  to  the  covariance  at  AR-0.  Again,  one  observes 
the  same  kind  of  contraction  in  the  covariance  as  noted  in  figure 
3.1.  [See  appendix  A. 2.] 

FIGURE  3.4 

This  is  the  wavenumber  analogue  of  the  frequency  spectnun  of 
figure  3.2,  now  with  and  is  obtained  from  (2.35a,b)  with 

A-/A„  ■  5**.  The  rectification  operation  similarly  spreads  the 
wavenumber  spectrum,  with  the  greatest  distortion  (v*0)  yielding 
the  greatest  wavenumber  spread,  as  expected  from  the 
corresponding  contraction  of  the  associated  covariance,  cf. 
figure  3.3  above.  [See  appendix  A. 4.] 
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CLASS  A  PLUS  GAUSS  NOISE 

FIGURE  3.5 

The  temporal  covariance  here  is  given  by  the  general  result 
(2.7),  with  the  associated  relations  (2.7),  (2.8),  (2.9),  wherein 
AR*0,  so  that  f  “T '■t2-tj^,  as  before,  and  where  (2.7a),  is 

now  given  analytically  by  (2.10)  for  v  «  0,1,2.  Here,  the 
parameter  values  are  *  AWg/p  “5^,  as  before,  now  with  A«0.2, 
r^*10“^,  AuL/p  s  A«j^  *  1  for  the  Class  A  non-Gaussian  noise 
component,  typically. 

Again,  for  the  super-clipper  (v“0),  the  contraction  in  the 
normalized  covariance  is  greatest,  cf.  figure  3.1.  But  the 
contribution  of  the  comparatively  strong  non-Gaussian  component 
exaggerates  this  effect.  [See  appendix  A.I.] 

FIGURE  3.6 

The  corresponding  intensity  (frequency)  spectrum  (AR=0), 
obtained  from  (2.7)  in  (2.39),  however,  shows  a  fine- structure 
not  exhibited  when  Gauss  noise  alone  (A»0)  is  applied  to  these 
ZMNL  devices.  The  spectral  levels  for  the  case  v-^O,  (A*0)  and 
(A>0),  cf.  figure  3.2  with  figure  3.6,  are  approximately  the 
same,  whereas  the  other  inputs,  cases  v»l,2,  are  much  elevated  as 
f  becomes  larger,  again  due  to  the  presence  of  the  structured 
Class  A  noise,  when  ^T^  i  1,  cf.  (2.3):  on  the  average,  the 
original  Class  A  "signals"  are  of  comparatively  short  duration, 
or  spectrally  wide  to  begin  with,  so  that  clipping  further 
spreads  the  spectrum.  [See  appendix  A. 3.] 

FIGURE  3.7 

The  spatial  covariance  when  Class  A  noise  is  added  to  the 
Gaussian  input  shows  analogous  behavior,  cf.  figures  3.3  and  3.5: 
the  covariance  is  compressed  vis-A-vis  the  input,  but  more  so 
than  in  the  Gauss-alone  situations.  Again,  (2.7)  -  (2.10)  are 
employed.  [See  appendix  A. 2.] 
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FIGURE  3.8 

The  corresponding  wavenumber  (intensity)  spectrum  with  Class 
A  noise  and  the  Gaussian  component,  obtained  from  (2.7)  -  (2.10) 
in  (2.39),  is  shown  here.  Comparison  with  figure  3.4  indicates  a 
broader  spectral  input,  due  to  the  non-Gaussian  component,  but  a 
relatively  narrower  output,  although  the  latter  is  still 
noticeably  spread  vis-A-vis  the  original  input.  [See  appendix 
A.4.  ] 

FIGURE  3.9 

Finally,  we  consider  the  angle-modulation  cases  described  in 

section  2.2-6  above,  where  weak  to  strong  angle  modulations 

-3 

(p  ~  1  to  50)  by  Class  A  noise,  with  a  weak  (F^^IO  )  Gaussian 
modulation  component,  is  employed. 

For  phase  modulation  by  non-Gaussian  noise,  based  on  (2.50) 
with  (2.44b),  (2.45),  (2.46),  the  resulting  normalized  intensity 
(frequency)  spectra  are  obtained  by  applying  (2.50)  to  (2.52), 
where  f  ■  u/2n;  u  ■  (w-Wq)/AWjj,  cf.  (2.49).  Note  the  "spike"  at 
f  ~  0.1,  followed  by  a  variety  of  sidelobes  which  rise  as  the 
phase  modulation  index  //p  increases.  The  spike  is  now  bounded  at 
f  “  0.8,  at  the  -10  dB  level,  when  Pp  *  50.  As  expected,  the 
larger  indexes  (Pp)  produce  broader  spectra.  [See  appendix  A. 5.] 

FIGURE  3.10 

For  frequency  modulation  by  non-Gaussian  noise,  from  (2.48) 
with  (2.49)  and  (2.44a),  the  corresponding  intensity  (frequency) 
spectra  again  exhibit  a  continuous  spike  (f  <  0.1).  With  small 
modulation  indexes  (Pp),  the  spectra  are  less  broad  than  for  the 
larger  indexes,  as  expected.  The  non-Gaussian  noise  component 
dominates  the  spectrum  here.  [See  appendix  A. 6.] 
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EXTENSIONS 


Other  situations  where  the  second-order  Class  A  probability 
density  functions  may  be  applied  are  noted  in  [2]  and  We 

list  some  of  the  extensions  of  the  analysis  to  the  following 
"classical"  problems: 

1)  The  inclusion  of  representative  signals,  with  Gauss  and 
non-Gauss  (Class  A)  noise,  in  the  problems  already  treated 
here  ( section  2 ) ; 

2)  The  case  of  the  full-wave  square-law  rectifier,  with  both 
Class  A  and  B  noise,  as  well  as  Gauss  noise; 

3)  The  extension  of  2)  to  include  general  broadband  and 
narrowband  signals; 

4)  The  calculation  of  signal-to-noise  ratios  and  deflection 
criteria,  cf.  [1;  section  5.3-4]. 

5)  Covariances  and  spectra  for  ZMNL  system  outputs,  with 
signals  as  well  as  non-Gaussian  noise  inputs; 

6)  The  rdle  of  the  electromagnetic  (or  acoustic)  interference 
(EMI  or  Acl)  scenario,  cf.  [5;  section  2B,5]; 

7)  Evaluation  of  the  large  (FM,PM)  indexes,  or  asymptotically 
Gaussian  cases,  cf.  [12]. 

Further  opportunities  to  extend  the  classical  theory  [2], [3], 
now  with  non-Gaussian  noise  inputs,  are  evident  from  the  examples 
and  methods  described  in  [1;  chapters  5,  12  -  16],  for  instance. 
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NORMALIZED  COVARIANCE 


NORMALIZED  FREQUENCY  f 


FIGURE  3.2  FREQUENCY  (INTENSITY)  SPECTRUM  (FOR  ^-0);  GAUSS 
NOISE  ONLY;  CF.  (2.39),  USED  WITH  (2.11),  AND  APPENDIX  A. 3 


NORMALIZED  SEPARATION  AR 


FIGURE  3.3  SPATIAL  COVARIANCE  (FOR  GAUSS 

NOISE  ONLY;  CF.  (2.11),  (2.7b),  AND  APPENDIX  A. 2 


SPECTRAL  LEVEL  (dB) 


NORMALIZED  COVARIANCE 
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FIGURE  3.5  TEMPORAL  COVARIANCE  (FOR  AR«0 ) ;  CLASS  A 
AND  GAUSS  NOISE;  CF.  (2.7)-(2.9)  AND  APPENDIX  A.l 
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NORMALIZED  SEPARATION  AR 


FIGURE  3.7  SPATIAL  COVARIANCE  (FOR  f',f-0);  CLASS  A 
AND  GAUSS  NOISE;  CF.  (2.7)-(2.10)  WITH  APPENDIX  A. 2 
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NORMALIZED  WAVENUMBER  Ic 

FIGURE  3.8  WAVENUMBER  SPECTRUM  (FOR  CLASS  A  AND 

GAUSS  NOISE;  CF.  (2.7)-(2.10)  IN  {2.38a,b)  AND  APPENDIX  A. 4 
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.01  .1  1  10  100 
NORMALIZED  FREQUENCY 

FIGURE  3.9a  PHASE  MODULATION  (INTENSITY)  SPECTRUM  FOR 
INDEX  /Jp«l,2,5,  CLASS  A  AND  GAUSS  NOISE;  CF.  (2.50)  WITH 
(2.44b),  (2.45),  (2.46)  IN  (2.52),  AND  APPENDIX  A. 5 
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FIGURE  3.9b  PHASE  MODULATION  (INTENSITY)  SPECTRUM  FOR 
INDEX  ^p»10,20,50,  CLASS  A  AND  GAUSS  NOISE;  CF.  (2.50)  WITH 
(2.44b),  (2.45),  (2.46)  IN  (2.52),  AND  APPENDIX  A. 5 
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.0001  .001  .01  .1  1  10  100 
NORMALIZED  FREQUENCY 

FIGURE  3.10a  FREQUENCY  MODULATION  (INTENSITY)  SPECTRUM 
FOR  INDEX  /jp-1,2,5,  CLASS  A  AND  GAUSS  NOISE; 

CF.  (2.48)  WITH  (2.49),  (2.44a),  AND  APPENDIX  A. 6 
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FIGURE  3.10b  FREQUENCY  MODULATION  (INTENSITY)  SPECTRUM 
FOR  INDEX  /ij.«10,20,50,  CLASS  A  AND  GAUSS  NOISE; 

CF.  (2.48)  WITH  (2.49),  (2.44a),  AND  APPENDIX  A. 6 


41/42 

Reverse  Blank 


TR  8887 


PART  II.  MATHEMATICAL  AND  COMPUTATIONAL  PROCEDURES 

4.  SOME  PROPERTIES  OF  THE  COVARIANCE  FUNCTION 

In  this  section,  we  collect  some  useful  relations  for  the 
covariance  and  auxiliary  functions  encountered  in  the  numerical 
evaluation.  These  are  necessary  for  rapid  computation  of  the 
multiple  series  involved  here  and  also  serve  as  checks  on  the 
numerical  procedures  employed. 

4.1  SIMPLIFICATION  AND  EVALUATION  OF  B^(Y) 

The  function  B^{y)  is  defined  by  the  following  combination  of 
hypergeometric  functions : 

B.(V)  -  f(-  I  -  i,  y")  * 

+  2  r2(^  .  l)  y  F(i-^,  y^)  for  y^*  i  1  .  (4.1) 

For  the  upper  F  function  in  (4.1),  we  have  [1;  (A. 1.39b)] 

V  «  0,  f(0,  0;  Y^)  -  1  ; 

V  -  1,  f(-  Y^)  -  Y  arcsin(Y)  +  (l  -  ; 

V  -  2,  f(-1,  -1;  Y^)  -  1  +  2Y^  ;  (4.2) 

where  arcs in  is  the  principal  value  inverse  sine  function.  On 
the  other  hand,  for  the  latter  F  function  in  (4.1),  we  have 
[1;  (A. 1.39a)  and  (A. 1.39c)] 
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V  =  2, 


V  =  0, 


V  = 


i.  3. 

2'  2' 


1,  f(o,  0; 

V^)  -  f(l  - 


Y^j  =  arosin(Y)  . 


2^^  1  +  2Y^ 

Y^J  •«•  -^y—  arcs3.n(Y). 


(4.3) 


When  these  quantities  are  substituted  in  the  above  expression 
for  B^(Y),  we  find  the  following  relatively  simple  relations: 

Bq(Y)  =  n  +  2  arcsin(Y)  , 

B^(Y)  =  Y  arcsin(Y)  +  (l  -  Y^]^  +  |y  , 

jj 

B2(Y)  =  (^  +  Y^)(|  +  arcsin(Y))  +  |y[i  -  Y^)  .  (4.4) 

These  three  quantities  can  be  computed  simultaneously  by  the 
following  very  compact  computer  coding  in  BASIC: 

Y2=Y*Y 

Sq=SQR(l.-Y2) 

T=ASN(Y)+1. 5707963267948966 

B0=T+T 

Bl=y*T+Sq 

B2=( .5+Y2)*T+1.5*Y*Sq  (4.5) 

Thus,  the  rather  formidable  expression,  above,  for  B^(Y)  can  be 
evaluated  by  the  use  of  just  one  square  root  and  one  arcs in  when 

V  =  0 ,  1 ,  2 . 

The  following  limiting  values,  which  are  obvious,  are  needed 
for  various  special  cases: 
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Bo(0) 

Bj{0) 

B2(0) 

B3(0) 

B4(0) 


1  , 
n/4  , 

1  , 

9n/16  , 


Bod) 

Bi(l) 

B2{1) 

B3(1) 

B4d) 


2n  , 

n  / 

3n/2  , 
15n/4  , 
105n/8 


(4.6) 


These  are  special  cases  of 


B^{0)  -  r 


B^(l)  =  2  n**  r(v  +  i)  , 


(4.7) 


(4.8) 


the  latter  following  from  [10;  (15.1.20)]. 

4.2  LIMITING  VALUES  OF  THE  COVARIANCE  FUNCTION 

The  covariance  function  at  normalized  separation  AR  and 
delay  f  is  given  by  (2.7)  as 


My(^/f)  *  exp[-A(2-p)] 


[A(l  -  P)] 


mi+m2 


mj^=0  ni2®0 


m  ^ * 


n«0 


n  +  m. 

V/ 

n  + 

^  1  p' 

A  ^  ^A 

A  +  ^A 

v/2 


B^(Y)  , 


(4.9) 


where 


p(f )  -  max{0,  1  -  |f I )  , 


(4.10) 


Y  -  Y(mj,m2,n)  -  ,  — 


X  k,  +  ri 

A  L  AG 


(n  +  m. 


^  ^k 


(4.11) 


45 


TR  8887 


kj^  ®  kj^(AR,f)  =  exp 


-  ar2  -  i 


(^“l1 

2  .2] 

p 

V  > 

T 

(4.12) 


Cq  *  kg{AR,f)  *  exp 


^12 
l^G 


1 

AWg' 

'  -2] 

2 

P 

T 

(4.13) 


The  functions  p,  kj^,  kg  can  be  replaced  by  other  functional 
dependencies,  if  desired.  The  function  B^(Y)  has  been 
considered  earlier  and  considerably  simplified  for  v  =  0,  1,  2. 

4.3  VALUE  AT  INFINITY 

As  AR  or  f  ±®,  then 

p  0,  kj^  0,  kg  -»  0,  Y  -»  0  .  (4.14) 

(If  |fl  remains  less  than  1  as  ^  tends  to  infinity,  then  p  does 

not  approach  zero;  this  nuance  has  been  discussed  elsewhere  in 
this  report.)  Then,  it  follows  that 


“y  "  E  H  simgr 

^  _ n  _ n  X  A 


mj»0  m2*'0 


f^l 

-4  +  ri 

V/  ^ 

’m~ 

+  r: 

A  A 

A  A 

.  4 

v/2 


B^(0) 


exp(-A)  ^  ^  (a  +  r^) 
m«0 


v/21 


(4.15) 


because  the  sum  on  n  can  be  terminated  with  the  n  «  0  term. 

The  sum  on  m  can  be  effected  in  closed  form,  for  v  »  0,  2,  4, 
etc.,  by  using  the  following  results: 


exp(A)  , 


(4.16) 
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^ A*"  — ■  a"* 

5u  '  Z_  oJr^-TTT  ' 


ni=0 


m=l 


m=0 


in»l 


V'  a”*  .  a”* 

/_  (m  -  2)  1  Z_  (m  -  l)i  ■ 
in=2  ni=l 


(A  +  A)  exp{A)  .  (4.18) 


There  follows 


My(«) 


-  tIi  * 

Bis  *  *  '•»)'] 


for  V  =  0 

for  V  =  2 

for  V  =  4 


(4.19) 


The  case  for  v  «  1  requires  a  numerical  summation,  once  A  and 
are  specified.  When  these  limiting  values  are  subtracted  from 
the  correlation  function,  we  obtain  the  covariance  function. 
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4.4  VALUE  AT  THE  ORIGIN 

For  AR  “  0,  f  =  0,  then 

P  *  1#  ^  '  (4.20) 

and 


My(0,0) 


exp(-A) 


00 


n»0 


(4.21) 


because  the  sums  on  m^^  and  m2  can  be  terminated  with  the  zero 
terms,  thereby  also  leading  to  Y  =  1. 

The  sum  on  n  can  be  accomplished  in  closed  form,  for 
V  *  0,  1,  2,  etc.,  by  using  results  given  earlier.  There  follows 


2n 


for  V  =  0) 


My(0,0) 


"(»  *  r') 


*  *  '■aI'I 


for  V  *=  1  >  . 

for  V  =  2 


(4.22) 
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PART  III.  APPENDICES  AND  PROGRAMS 

APPENDIX  A.l  -  EVALUATION  OF  COVARIANCE  FUNCTION 
FOR  ZERO  SEPARATION  -  0) 

A  program  for  the  numerical  evaluation  of  covariance 
M  (AR,t)  for  AR  «  0  IS  contained  in  this  appendix.  Inputs 

y  2  2 

required  of  the  user  are  A,  ,  (A«g/p)  ,  &{t),  N(t)/ 

in  lines  20  -  70.  Since  we  are  generally  interested  in  values  of 

A  less  than  1,  the  series  for  in  (4.9)  will  not  have  to  be 

taken  to  very  large  values  of  m, ,  m~,  n;  accordingly,  the  values 
k  1  ^ 

of  {A  /kl )  are  tabulaLed  once  in  lines  260  -  300  with  a  tolerance 
of  lE-10  set  in  line  80. 

The  values  of  the  covariance  at  infinity,  as  given  by  (4.19), 
are  computed  and  subtracted  in  lines  220  -  240  and  400  -  420; 
this  is  in  anticipation  of  taking  a  Fourier  transform  of  a 
covariance  function  which  decays  to  zero  for  large  arguments  AR. 

The  functions  B^(Y)  and  Ay(AR,f )  are  available  in  the  two 
subroutines  starting  at  lines  1010  and  1120,  respectively.  The 
latter  subroutine  actually  calculates  the  covariance  at  general 
nonzero  values  of  both  AR  and  f,  although  we  only  employ  it  for 
AR  «  0  in  this  appendix;  see  lines  10  and  380.  Also,  for  AR  «  0, 
the  parameter  Lg2  *  (A^/A^)  is  not  relevant  and,  hence,  is 
entered  as  zero  in  line  380. 

The  exponential  Gaussian  forms  for  k^  and  k^  are  used  in 
lines  1200  and  1210,  while  the  triangular  form  for  p  is  entered 
in  line  1240.  Any  of  these  can  be  replaced,  if  desired,  by  forms 
more  appropriate  to  the  user. 

The  program  is  written  in  BASIC  for  the  Hewlett  Packard  9000 
Computer  Model  520.  The  designation  DOUBLE  denotes  integer 
variables,  not  double  precision.  The  output  from  the  program  is 
stored  in  data  files  AOTO,  AOTl,  AOT2,  for  v  ■  0,  1,  2, 
respectively. 
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10 

Rc=0. 

! 

DelR-' 

20 

fl=.  2 

! 

R<subA) 

30 

Gp^.001 

1 

GAMMA" (subfl) 

40 

N1b2=l . 

j 

(Del  UCsubD/Bet 

a)"'2 

50 

Wgb2=25. 

1 

(De1U(subG)/Bet 

a)^2 

60 

Dtc=.01 

1 

INCREMENT  IN  Tau-' 

70 

Ntc=200 

! 

NUMBER  OF  Tau'' 

VALUES 

80 

T  o 1 eranc  e= 1 . E  - 1 0 

90 

COM  fif<0;40),C<0!80),Sq< 

e:80> 

100 

COM  DOUBLE  J 

! 

INTEGER 

110 

DIM  Kag<200) , Tc  <0; 200> , F0<0: 200) , FI  (0! 200>  ,  F2 

<0:200 

120 

DOUBLE  Nlc,K 

1 

INTEGERS 

130 

FOR  K=0  TO  Htc 

140 

Tc=K*Dtc 

1 

Tau"' 

150 

Rho=MflX<0. , l.-fiBS<Tc)) 

1 

Rho 

160 

T2=.5*Tc*Tc 

170 

Kl=EXP<-Wlb2»T2) 

180 

Kg=EXP<-Wgb2*T2) 

190 

Kag<K)  =  <Rho*K1  +Gp*Kg)/<:  1 

•  +Gp> 

!  INPUT  COVARIANCE 

200 

NEXT  K 

210 

fll=l . /fl 

! 

A>0  REQUIRED 

220 

F0inf=PI 

230 

FI inf=FNFl inf <fl,Gp) 

240 

F2inf  =  .  25*PI*<  1 .  ■»-Gp)»<  1 . 

+Gp) 

250 

flf <0)=1. 

260 

FOR  K=1  TO  40 

270 

J=K 

280 

flf <K>=T*flf <K-1 )*fl/K 

1 

A'^K/K! 

290 

IF  T<Tolerance  THEN  320 

300 

NEXT  K 

310 

PRINT  ''40  TERMS  IN  flf<*) 

M 

320 

FOR  K=0  TO  J*2 

330 

C<K)=T=K*fll+Gp 

340 

Sq<K)=l./SQR<T) 

350 

NEXT  K 

360 

FOR  K=0  TO  Ntc 

370 

Tc<K)=Tc=K»Dtc 

! 

Tau^ 

380 

CALL  M>;c<Rc,Tc,fl,Gp,  Wlb2 

,Ugb2 

,0.  ,F0<K),F1<K), 

F2<K)) 

390 

NEXT  K 

400 

MAT  F0*F0-<F01nf ) 

410 

MAT  F1=F1-<F1 inf ) 

420 

MAT  F2=F2-<F2inf ) 

430 

MAT  F0=F0/<F0<0)) 

440 

MAT  F1=F1/<FI<0)) 

450 

MAT  F2*=F2/<F2<0)) 

460 

PRINT  “INFINITYj ";F0inf 5 

FI  inf 

}F2inf 

470 

PRINT  "MINIMA;  ";MIN<F0 

<♦)); 

MIN<F1<*>);MIN<F2<*)) 

480 

PRINT  "AT  Ntc;  ";F0<Ntc 

);F1< 

Ntc>}F2(Ntc) 

490 

CREATE  DATA  "flITl",8 

500 

ASSIGN  ttl  TO  "AITl" 

510 

PRINT  #ljKag<*> 

520 

CREATE  DATA  "AOT0",8 

530 

ASSIGN  #1  TO  "AOT0" 

540 

PRINT  <»ljF0<*> 

550 

CREATE  DATA  "AOTr‘,8 
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560  ASSIGN  #1  TO  “flOTl" 

570  PRINT 

580  CREATE  DATA  •■A0T2",8 

590  ASSIGN  #1  TO  "A0T2‘' 

600  PRINT  #ljF2(*) 

610  ASSIGN  ttl  TO  * 

620  Tcmax=Dlc*Ntc 

630  GINIT  200/260 

640  PLOTTER  IS  505, "HPGL" 

650  PRINTER  IS  505 

660  LIMIT  PLOTTER  505,0,200,0,260 

670  VIEWPORT  22,85,19,122 

680  WINDOW  0. , 1. ,0. , 1. 

690  PRINT  “735" 

700  GRID  .25,. 25 

710  PRINT  "VS36“ 

720  PLOT  Tc<»),Kag<*) 

730  PENUP 

740  PLOT  Tc<*),F0<*) 

750  PENUP 

760  PLOT  Tc<*),Fl<*) 

770  PENUP 

780  PLOT  Tc<*),F2(*) 

790  PENUP 

800  PAUSE 

810  PRINTER  IS  CRT 

820  PLOTTER  505  IS  TERMINATED 

830  END 

840  ! 

850  DEF  FNFlinf <R,Gp)  !  for  w<=nu)  =  1 

860  Tol=l.E-18 

870  Ag=A*Gp 

880  T=l. 

890  S*SQR<1.+Ag) 

900  FOR  M=2  TO  100 

910  T*T*A/M 

920  P«T*SQR<M+Ag> 

930  S=S+P 

940  IF  P<S»Tol  THEN  970 

950  NEXT  M 

960  PRINT  "100  TERMS  IN  FNFlinf" 

970  T=Gp+H*S»S+2.*SQR<Hg)*S 

980  RETURN  EXP<-2.*A)*T 

990  FNEND 

1000  ! 

1010  SUB  Bnu<Y,B0,Bl,B2)  !  Bv<Y>  for  v=0,l,2 

1020  IF  Y>1.  THEN  PRINT  "Y  =  1  +";Y-1. 

1030  IF  Y>1.  THEN  Y**!. 

1040  Y2»Y*Y 

1050  Sq=SQR<l.-Y2) 

1060  T*=ASN(Y)  +  1. 5707963267948966 

1070  B0>=T  +  T 

1080  B1>»Y*T  +  Sq 

1090  B2=<.5+Y2)*T+1.5*Y*Sq 

1100  SUBEND 

1110  I 
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1120  SUE  MMc<Rc,Tc,ft,Gp,Wlb2,Wgb2,Lg2,S0,Sl,S2> 

1130  COM  flf , C(*) , Sq(*> 

1140  COM  DOUBLE  J  !  INTEGER 

1150  ftLLOCflTE  flp<0: J) , flpl <0; J) 

1160  DOUBLE  K,M1,M2,N,K1,K2  !  INTEGERS 

1170  ftl=l./fl  !  fi>0  REQUIRED 

1180  T2=.5»Tc*Tc 

1190  R2=Rc»Rc 

1200  K1=EXP<-R2-Wlb2*T2) 

1210  Kg=EXP<-Lg2*R2-Wgb2*T2) 

1220  ftk=fll*Kl 

1230  Gk=Gp*Kg 

1240  Rho=NflX<0. , 1 . -flBSCTc ) )  !  Rho 

1250  Rho 1=1. -Rho 

1260  fip<0)=flpl<0)=Pk=Pkl=l. 

1270  FOR  K=1  TO  J 

1280  Pk=Pk*Rho 

1290  Pkl=Pkl*Rhol 

1300  T=flf<K> 

1310  flp<K)=T*Pk 

1320  flpl<K>=T»Pkl 

1330  NEXT  K 

1340  S0ml=Slml=S2ml=Q. 

1350  FOR  M1=0  TO  J 

1360  S0m2=Slni2=S2in2=0. 

1370  FOR  M2=0  TO  J 

1380  S0n=Sln»S2n=0. 

1390  FOR  N«0  TO  J 

1400  K1=N+M1 

1410  K2=N+M2 

1420  T=fip(N) 

1430  P=C<K1 )*C<K2) 

1440  Y=<N*flk+Gk )*Sq<Kl )*Sq<K2) 

1450  CALL  Bnu<Y, B0, B1 , B2) 

1460  S0n=S0n+T*B0 

1470  Sln=Sln+T*SQR<P)*Bl 

1480  S2n=S2n+T*P*B2 

1490  NEXT  N 

1500  T2=flpl<M2) 

1510  S0m2  =  S0iri2+T2*S0n 

1520  Sl(n2  =  Slm2  +  T2*Sln 

1530  S2ni2  =  S2iri2+T2*S2n 

1540  NEXT  M2 

1550  Tl=flpl<Ml) 

1560  S0ml=S0ml+Tl»S0m2 

1570  Slrrtl=Slml+Tl*Sltti2 

1580  S2nil=S2(nl+Tl*S2M2 

1590  NEXT  Ml 

1600  T=EXP<-fi*<2. -Rho) ) 

1610  S0=T*S0ml 

1620  Sl=T»Slml 

1630  S2  =  T*S2titil 

1640  SUBEND 
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APPENDIX  A. 2  -  EVALUATION  OF  COVARIANCE  FUNCTION 
FOR  ZERO  DELAY  (f,f'  =  0) 


A  program  for  the  numerical  evaluation  of  covariance 
M  (AR,f )  for  f,f'  «  0  is  contained  in  this  appendix.  Inputs 

y  2  ^  ^ 

required  of  the  user  are  A,  ,  &(AR),  N(AR),  in  lines 

20  -  60.  The  tolerance  for  terminating  the  triple  infinite  sums 
is  set  at  IE- 15  in  line  70.  The  output  from  the  program  is 
stored  in  data  files  AORO,  AORl,  AOR2,  for  v  «=  0,  1,  2, 
respectively.  Other  relevant  comments  are  made  in  appendix  A.l. 

The  limit  of  at  AR  =  ®  (when  f  =0)  is  given  by  the  closed 
form  results 


n 


for  V  =  0 


My(“,0) 


1  +  r; 

A 


Ilk  * 


for  V  =  1 

for  V  =  2 


(A. 2-1) 


A 

These  values  have  been  subtracted  from  M  so  that  we  can  Bessel 

y  A* 

transform  a  function  which  tends  to  zero  as  AR 


le 

20 

30 

40 

50 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 


Tc=0. 
fl=.  2 
Gp=. 001 
Lg2=5. 

Drc=. 005 
Nrc  =900 

To  I erance=l . E-15 


!  Tau'' 

!  fiCsubfl) 

!  GfiMMfl' (subR> 

!  <DelL/De1G)''2 
!  INCREMENT  IN  DclR-^ 

!  NUMBER  OF  DelR-  VALUES 


COM  flf<0:40),C<0;80),Sq<0;80) 

COM  DOUBLE  J  !  INTEGER 

DIM  Rc  <0; 900) , Kag<0: 900) , F0<0: 900), FI <0: 900) , F2<0: 900) 
DOUBLE  Nrc,K  !  INTEGERS 

Al=l./R  !  R>0  REQUIRED 

F0inf=PI  !  LIMITS  FOR 

Flinf=l.+Gp  !  Rc=infinity 

F2irif=.25»PI*(<l.+Gp)*<l.+Gp>  +  Rl)  !  RND  Tc=0 
Rf (0)=1 . 

FOR  K=1  TO  40 
J=K 


flf  <K)=T  =  Rf  <K-1  )*H/'K  !  R^'K/'K! 

IF  T<Tolerance  THEN  230 
NEXT  K 

PRINT  "40  TERMS  IN  RfC*)“ 

FOR  K«0  TO  J*2 
C<K)=T«K*R1+Gp 
SqCK)«l . /SORCT) 

NEXT  K 
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270  FOR  K=0  TO  Nrc 

280  Rc  <K)=Rc=K*Dt'C  !  DelR-' 

290  R2=Rc«Rc 

300  K1=EXP<-R2> 

310  Kg=EXP<-Lg2*R2) 

320  Rho=MaX(0. , 1 . -fiBS<Tc ) )  !  Rho 

330  Kag<K)=<Rho*Kl+Gp*Kg)/<l.+Gp>  !  INPUT  COVARIANCE 

340  CALL  rivc<Rc,Tc,A,Gp,0.  ,0.  ,Lg2,F0<K>,Fl(K),F2(K>) 

350  NEXT  K 

360  MAT  F0=F0-<F0inf ) 

370  MAT  Fl=Fl-<:Flinf) 

380  MAT  F2=F2-<F2inf ) 

390  MAT  F0=F0/'<F0<0>  ) 

400  MAT  F1=F1/<F1<0)> 

410  MAT  F2  =  F2/'<F2<0)) 

420  PRINT  “INFINITY: “;F0 inf} FI  inf; F2 inf 

430  PRINT  "MINIMA:  " ; MIN<F0<*) > ; MIN<F1 <*> ) ; MIN<F2<*) > 

440  PRINT  "AT  Nrc:  "  ; F0 < Nrc ) } F 1 ( Nrc > ; F2 < Nrc > 

450  CREATE  DATA  "AIR1",33 

460  ASSIGN  #1  TO  "AIRl" 

470  PRINT  #l;Kag<*) 

480  CREATE  DATA  "AOR0",33 

490  ASSIGN  #1  TO  "AOR0" 

500  PRINT  #1;F0<*) 

510  CREATE  DATA  "A0R1",33 

520  ASSIGN  #1  TO  "AORl" 

530  PRINT  #l}Fl<*i 

540  CREATE  DATA  "R0R2",33 

550  ASSIGN  #1  TO  "R0R2" 

560  PRINT  #1;F2<*) 

570  ASSIGN  #1  TO  * 

580  Rcmax=Drc*Nrc 

590  GINIT  200/260 

600  PLOTTER  IS  505, "HPGL" 

610  PRINTER  IS  505 

620  LIMIT  PLOTTER  505,0,200,0,260 

630  VIEWPORT  22,85,19,122 

640  WINDOW  0. ,3. ,0. , 1. 

650  PRINT  "VS5" 

660  GRID  .5, .25 

670  PRINT  "VS36" 

680  PLOT  Rc<*),Kag<*) 

690  PENUP 

700  PLOT  Rc<*),F0<*) 

710  PENUP 

720  PLOT  Rc<*),Fl<*> 

730  PENUP 

740  PLOT  Rc(*),F2<*) 

750  PENUP 

760  PAUSE 

770  PRINTER  IS  CRT 

780  PLOTTER  505  IS  TERMINATED 

790  END 

800  I 

810  SUB  Bnu<Y,B0,Bl ,B2>  !  Bw<Y)  for  o=0,l,2 

820  !  SEE  APPENDIX  A.  1 

900  SUBEND 

910  ! 

920  SUB  Myc<Rc,Tc,A,Gp,W1b2,Mgb2,Lg2,S0,Sl,S2> 

930  I  SEE  APPENDIX  A. 1 

1440  SUBEND 
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APPENDIX  A. 3  -  EVALUATION  OF  TEMPORAL  INTENSITY 

SPECTRUM  FOR  ZERO  SEPARATION  (AR  =  0) 

A  program  for  the  numerical  evaluation  of  the  Fourier 
transform  of  covariance  My(0,f)  -  My(®)  is  contained  in  this 
appendix.  Inputs  required  of  the  user  are  listed  in  lines 
10  -  30.  The  data  input,  AOTO  or  AOTl  or  AOT2,  as  generated  by 
means  of  the  program  in  appendix  A.l,  is  injected  by  means  of 
lines  410,  600,  and  790. 

In  order  to  keep  the  FFT  (fast  Fourier  transform)  size,  N  in 
lines  30  and  320,  at  reasonable  values,  the  data  sequence  is 
collapsed,  without  any  loss  of  accuracy,  according  to  the  method 
given  in  [8;  pages  7-8]  and  [9;  pages  13  -  16].  The 
integration  rule  documented  here  is  the  trapezoidal  rule;  this 
procedure  is  very  accurate  and  efficient  and  is  recommended  for 
numerical  Fourier  transforms. 


10 

Ntc=200 

!  NUMBER  OF 

Tau'^  VALUES 

20 

Dtc=. 01 

1  INCREMENT 

IN  Tau'^ 

30 

H=1024 

!  SIZE  OF  FFT;  \i  >  Ntc  REQUIRED 

40 

DOUBLE  Ntc,N,N4,H2,Ns 

!  INTEGERS 

50 

N4=N/4 

60 

N2=N/'2 

70 

REDIM  Cos<0;N4),X<0;N-l) 

,Y<0;N-1) 

80 

DIM  Cos<256) ,X< 1023) , Y< 1 023 ) , fl<200 ) 

90 

T=2. »PI/N 

100 

FOR  Ns=0  TO  N4 

110 

Cos<Ns)=COS<T*Ns) 

!  QUARTER-COSINE  TABLE  IN  Cose*) 

120 

NEXT  Ns 

130 

GINIT  2O0/'260 

140 

PLOTTER  IS  505, "HPGL" 

150 

PRINTER  IS  505 

160 

LIMIT  PLOTTER  505,0,200, 

0,260 

170 

VIEWPORT  22,85, 19, 122 

180 

WINDOW  0,N2,-5, 1 

190 

PRINT  "VS5" 

200 

GRID  N^ie, 1 

210 

PRINT  "VS36" 

220 

fiSSIGN  ttl  TO  "fllTl" 

230 

REND  ttl;n(*) 

240 

MAT  X=<0. > 

250 

HAT  Y=<0. ) 

260 

X<0)=. 5*A<0) 

270 

FOR  Ns=l  TO  Ntc-1 

280 

X<Ns)«A<Ns) 

290 

NEXT  Ns 

300 

X<Nlc)*.5*A<Ntc> 
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310  MftT  X  =  X*(Dt.c*4.  ) 

320  CALL  Fft  14-:N,Co£<*),X<:*),Y<*)) 

330  FOR  H£=0  to  N2 

340  fir=X<Ns) 

350  IF  ftr>0.  THEN  380 

360  PENUP 

370  GOTO  390 

380  PLOT  Ns,LGT<flr) 

390  NEXT  Ns 

400  PENUP 

410  ASSIGN  «1  TO  "AOT0" 

420  READ  1I1}A<*) 

430  MAT  X=<0.> 

440  MAT  Y=<0.) 

450  X<0)=.5*R(0) 

460  FOR  Ns=l  TO  Ntc-1 

470  X(N£)=A<Ns) 

480  NEXT  Ns 

490  X<Ntc )=. 5*A<Ntc ) 

500  MAT  X  =  X»(Dtc*4.  ) 

510  CALL  Fft 14<H,Co£<*),X<*>,Y<*)) 

520  FOR  Ns=0  TO  N2 

530  Ar=X<Ns) 

540  IF  Ar>0.  THEN  570 

550  PENUP 

560  GOTO  580 

570  PLOT  Ns,LGT<Ar) 

580  NEXT  Ns 

590  PENUP 

600  ASSIGN  #1  TO  ''AOTl" 

610  READ  ttliA<*> 

620  NAT  X=<0.) 

630  MAT  Y=<0. ) 

640  XC0)=.5*A<0) 

650  FOR  Ns=l  TO  Ntc-1 

660  X(Ns)=A<Ns) 

670  NEXT  Ns 

680  X<Ntc )«. 5*A<Ntc ) 

690  MAT  X*X*<Dtc*4. ) 

700  CALL  Fftl4<N,Cos<»),X<*),Y<*)) 

710  FOR  Ns«0  TO  N2 

720  Ar=X<Ns> 

730  IF  Ar>0.  THEN  760 

740  PENUP 

750  GOTO  770 

760  PLOT  Ns,LGT<Ar) 

770  NEXT  Ns 

780  PENUP 

790  ASSIGN  ttl  TO  "A0T2” 

800  READ  ttl|A<*> 

810  MAT  A=A/<A<0)) 

820  MAT  X=<0. ) 

830  MAT  Y=<0.> 
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840  X<0)=.5*fl<0) 

850  FOR  Ns=l  TO  Ntc-1 

860  X<H£)=fl<Ns) 

870  NEXT  Ns 

880  X<Ntc )=. 5*fi<Ntc > 

890  MAT  X=X*(Dtc*4. > 

900  CALL  Fft 14<N,Cos<*),X<*),Y<*)) 

910  FOR  Ns=0  TO  N2 

920  Ar=X<Ns> 

930  IF  Ar>0.  THEN  960 

940  PENUP 

950  GOTO  970 

960  PLOT  Ns,LGT<Ar) 

970  NEXT  Ns 

980  PENUP 

990  PAUSE 

1000  END 

1010  ! 

1020  SUE  Fftl4<D0UBLE  N.REAL  Cos<*)  ,  X<»>  ,  Y<*>  >  !  N<  =2'' 1  4=  1 6384  ;  0  SUES 

1030  DOUBLE  Log2n,Nl,N2,N3,N4, J,K  !  INTEGERS  <  2^31  =  2,147,483,648 

1040  DOUBLE  II, 12, 13, 14, 15, 16, 17, 18, 19, 110, Ill, 112, 113, 114, Lc:0:13> 

1050  IF  N=1  THEN  SUBEXIT 

1060  IF  N>2  THEN  1140 

1070  A=X<0)+X<1) 

1080  X< 1 )=X<0)-X<  1) 

1090  X(0)=A 

1100  R=Y<0)+Y<n 

1110  Y<n=Y<0)-Y<l) 

1120  Y<0)=A 

1130  SUBEXIT 

1140  A=L0G<N)/L0G<2.  ) 

1150  Log2n=A 

1160  IF  HBS<A-Log2n)<l.E-8  THEN  1190 

1170  PRINT  “N  ='‘jN;‘'IS  NOT  A  POWER  OF  2;  DISALLOWED." 

1180  PAUSE 

1190  Nl=N/'4 

1200  N2=N1+1 

1210  N3=N2+1 

1220  N4=N3+N1 

1230  FOR  11=1  TO  Log2n 

1240  I2  =  2^<Log2n-1 1  > 

1250  I3=2«I2 

1260  I4=N/I3 

1270  FOR  15=1  TO  12 

1280  I6=( 15-1 )*I4+1 

1290  IF  I6<=N2  THEN  1330 

1300  Rl=-Cos<N4-I6-l ) 

1310  r;2  =  - Cos<16-Nl-l) 

1320  GOTO  1350 

1330  Rl=Cos<I6-l) 

1340  A2=-Co£<N3-I6-1 ) 

1350  FOR  17=0  TO  N-I3  STEP  13 

1360  18=17+15-1 

1370  19=18+12 

1380  T1=XCI8) 

1390  T2=X<I9) 
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1400 

T3=YCI8> 

1416 

T4=Y( I9> 

1420 

fl3=Tl-T2 

1430 

fl4=T3-T4 

1440 

X< I8)=T1+T2 

1450 

Y< I8)=T3+T4 

1460 

X< I9)=fll*fl3-fl2*fl4 

1470 

Y< I9)=fll*fl4+fl2*fl3 

1480 

NEXT  17 

1490 

NEXT  15 

1500 

NEXT  11 

1510 

1 1 =Log2n+ 1 

1520 

FOR  12=1  TO  14 

1530 

L< 12-1 )=1 

1540 

IF  I2>Log2n  THEN  1560 

1550 

L<  12-n=2'<  11-12) 

1560 

NEXT  12 

1570 

K  =  0 

1580 

FOR  11=1  TO  L<13) 

1590 

FOR  12=11  TO  L<12) 

STEP  L<13) 

1600 

FOR  13=12  TO  L<11) 

STEP  L<12) 

1610 

FOR  14=13  TO  L<10) 

STEP  L<11) 

1620 

FOR  15=14  TO  L(9)  STEP  L<10) 

1630 

FOR  16=15  TO  L<8)  STEP  L<9) 

1640 

FOR  17=16  TO  L<7)  STEP  L<8) 

1650 

FOR  18=17  TO  L<6)  S 

TEP  L<7) 

1660 

FOR  19=18  TO  L<5)  STEP  L<6) 

1670 

FOR  110=19  TO  L<4) 

STEP  L<5) 

1680 

FOR  111=110  TO  L<3) 

STEP  L<4> 

1690 

FOR  112=111  TO  L<2) 

STEP  L<3) 

1700 

FOR  113=112  TO  L<1) 

STEP  L<2) 

1710 

FOR  114=113  TO  L<0) 

STEP  L<1) 

1720 

J=I 14-1 

1730 

IF  K>J  THEN  1800 

1740 

fi=X<K) 

1750 

X<K)=X< J) 

1760 

X< J)=fl 

1770 

fl  =  Y<K) 

1780 

Y<K)=Y<:J) 

1790 

Y< J)=fl 

1600 

K=K+1 

1810 

NEXT  114 

1820 

NEXT  113 

1830 

NEXT  112 

1840 

NEXT  Ill 

1850 

NEXT  110 

1860 

NEXT  ;9 

1870 

NEXT  18 

1880 

NEXT  17 

1890 

NEXT  16 

1900 

NEXT  15 

1910 

NEXT  14 

1920 

NEXT  13 

1930 

NEXT  12 

1946 

NEXT  I  1 

1950 

SUBENO 
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APPENDIX  A. 4  -  EVALUATION  OF  WAVENUMBER  INTENSITY 
SPECTRUM  FOR  ZERO  DELAY  (f,f'  =  0) 


A  program  for  the  numerical  evaluation  of  the  zeroth-order 
Bessel  transform  of  covariance  My(AR,0)  -  My(<»>)  is  contained  in 
this  appendix.  Inputs  required  of  the  user  are  listed  in  lines 
10  -  40  and  are  coupled  to  appendix  A. 2,  where  the  data  input, 
AORO  or  AORl  or  AOR2 ,  was  generated.  The  numerical  Bessel 
transform  is  accomplished  by  means  of  Simpson's  rule  with  end 
correction  [11;  pages  414  -  418],  and  is  exceedingly  accurate  for 
the  small  increment,  .005,  in  AR  employed  in  line  30. 


10 

20 

30 

40 

50 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 


Dkc=. 4 
Hkc=200 
Drc  = . 005 
Hrc  =900 

DOUBLE  Hrc,Nkc,I,H£ 

REDIM  C<0:Hrc) 

REDIM  Hi <0;Nkc),U0<0;Hk 
DIM  C<900) , Hi <200) , W0<200 
ASSIGN  #1  TO  "flIRl" 

READ  #1;C<>) 

FOR  1=0  TO  Nkc 

Kc=I*Dkc  ! 

T=Kc*Drc 

Se=So=0 . 

FOR  Ns=l  TO  Nrc-l  STEP  2 
So=So+N£*FNJo<T*M£)*C<Hs) 
NEXT  Ns 

FOR  N£=2  to  Nrc-2  STEP  2 
Se=Se+Ns*FNJo<T*N£)*C<Ns) 
NEXT  Ns 

Hi < I  )=C<0)  +  16. «So+14. *Se 
NEXT  I 

MAT  Hi  =Wi  »<Drc*Drc*2.  *PI/' 
ASSIGN  #1  TO  "AORO" 

READ  #ljC<*) 

FOR  1=0  TO  Nkc 
Kc  =  I *Dk  c 
T=Kc*Drc 
Se=So=0 • 

FOR  Ns=l  TO  Nrc-l  STEP  2 
So=So+N£*FNJo<T*Ns)*C<Ns) 
NEXT  Ns 

FOR  Ns=2  TO  Nrc-2  STEP  2 
Se=Se+N£*FNJo<T*Ns)*C<N£) 
NEXT  Ns 


INCREMENT  IN  k- 
NUMBER  OF  k^  VALUES 
INCREMENT  IN  DelR^ 
NUMBER  OF  DelR'^  VALUES 
INTEGERS 

Wl<0:Nkc),W2(0:Nkc> 

,W1<200),W2<200> 


W0< I )=C<0>+16. *80+14. *Se 
NEXT  I 

MAT  W0  =  W0*<Drc*Drc*2.  *PI/'15.  ) 
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390  ASSIGN  #1  TO  "RORl" 

400  REAL  #l5Ci;*) 

410  FOR  1=0  TO  Nkc 

420  Kc=I*Dkc 

430  T=Kc*Drc 

440  Se=So=0. 

450  FOR  Hs=l  TO  Nrc-1  STEP  2 

460  So=So+N£*FNJo<T*Ns)*C<Hs> 

470  NEXT  Ns 

480  FOR  Ns=2  TO  Nrc-2  STEP  2 

490  Se=Se+Ns*FNJo<T*N£)*C<Ns) 

500  NEXT  Ns 

510  W1 < I )=C<0)+16. *So+14. *Se 

520  NEXT  I 

530  MAT  Wl=Wl*<Drc*Drc*2.  *PI/'15.  > 

540  ASSIGN  #1  TO  ■•A0R2" 

550  READ  #1;C<*) 

560  ASSIGN  #1  TO  * 

570  FOR  1=0  TO  Nkc 

580  Kc=I*Dkc 

590  T=Kc*Drc 

600  Se=So=0i 

610  FOR  N£=1  to  Nrc-1  STEP  2 

620  So=So+Ns*FNJo<T*Ns)*C<Ns> 

630  NEXT  Ns 

640  FOR  Ns=2  TO  Nrc-2  STEP  2 

650  Se=Se+Ns*FNJo<T«Ns)*C<Ns) 

660  NEXT  Ns 

670  W2< I >=C<0)+16. *So+14.*Se 

680  NEXT  I 

690  MAT  W2=W2*<Drc*Drc*2. *PI/15. ) 

700  GINIT  200/260 

710  PLOTTER  IS  505, "HPCL" 

720  PRINTER  IS  505 

730  LIMIT  PLOTTER  505,0,200,0,260 

740  VIEWPORT  22,85,19,122 

750  WINDOW  0, Nkc, -9, 1 

760  PRINT  "VSS" 

770  GRID  25, 1 

780  PRINT  "VS36" 

790  FOR  1=0  TO  Nkc 

800  W=Wi(I> 

810  IF  W>0.  THEN  840 

820  PENUP 

830  GOTO  850 

840  PLOT  I,LGT<W) 

850  NEXT  I 

860  PENUP 
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870 

FOR  1=0  TO  Nkc 

880 

N=N0< I > 

890 

IF  W>0.  THEN  920 

900 

PENUP 

910 

GOTO  930 

920 

PLOT  I,LGT<W) 

930 

NEXT  I 

940 

PENUP 

950 

FOR  1=0  TO  Nkc 

960 

W=W1 <  I  ) 

970 

IF  W>0.  THEN  1000 

980 

PENUP 

990 

GOTO  1010 

1000 

PLOT  I,LGT<W) 

1010 

NEXT  I 

1020 

PENUP 

1030 

FOR  1=0  TO  Nkc 

1040 

W=W2<  I  ) 

1050 

IF  W>0.  THEN  1080 

1060 

PENUP 

1070 

GOTO  1090 

1080 

PLOT  I,LGT<W) 

1090 

NEXT  I 

1  100 

PENUP 

1 110 

PAUSE 

1120 

PRINTER  IS  CRT 

1130 

PLOTTER  505  IS  TERMINATED 

1140 

END 

1150 

! 

1160 

DEF  FNJo<X)  !  Jo<X) 

FOR  ALL  X 

1170 

Y=ABS<X) 

1180 

IF  Y>8.  THEN  1280 

1190 

T=Y*Y  !  HART, 

#5845 

1200  P  =  2271490439. 553€033-T*< 551  3584. 5647707522-1*5292. 6171303845574  ;' 

1210  P  =  2334489171877869. 7-T*<47765559442673. 588-T* (462 172225031 . 7 1803-T*F  >  > 
1220  P=1 8596231 762 1897804. -T*< 441 45829391 815982. -T*P) 

1230  Q  =  20425 1483. 52 134357  +  T*<494030. 79491 81 3972  +  T* <884. 72036756 175504  +  T>  > 

1240  0=2344750013658996. 8+T*< 15015462449769. 752  +  T* < 64398674535 . 133256  +  T*Q>  > 

1250  0=185962317621897733. +T»Q 

1260  Jo=P/0  ' 

1270  RETURN  Jo 

1280  Z  =  8.-'Y  !  HART,  #6546  8.  6946 

1290  T=2*2 

1300  Pn  =  2204. 5010439651804  +  T*< 128. 67758574871 4 19+T*. 9004  793474S02S803> 

1310  Pn=8554. 82254 150666 17+T*< 8894. 4375329606 194+T*Pn) 

1320  Pd=2214.0488519147104+T»< 130.88490049992388+T) 

1330  Pd=8554. 82254 15066628  +  T*< 8903. 836 1 4 1 7095954  +  T*Pd > 

1340  0n=13. 990976865960680+T»<l .0497327982345548+T*. 0093525953294O319> 

1350  0n=-37. 51 0534954957 1  12-T* <46. 0938268 14625 175+T*0ri) 

1360  0d=921 .56697552653090+T*<74. 42838974 141 1 179+T) 

1370  Qd=2400. 6742371 172675+T*< 2971 . 9837452084920+T*Qd > 

1380  T=Y-. 78539816339744828 

1390  Jo=.  232094791 77387820*SQR<Z>*<COS<T>*Pn/'Pd-S I H<T>*2*Qri/C!d> 

1400  RETURN  Jo 

1410  FNEND 
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APPENDIX  A. 5  -  EVALUATION  OF  PHASE  MODULATION  INTENSITY  SPECTRUM 

The  normalized  covariance  function  for  phase  modulation  is 
given  by  (2.50)  in  the  main  text,  namely 

»  expf-  r;  fjl  fl  -  exp(-C)]  -  A[2  -  p(C)]  + 

^  ^  ^  (A. 5-1) 

^2 

+  2A  [1  -  p(C)]  exp[-p2/A)  +  A  p(C)  exp(-  [1  -  exp(-bC)])] 

for  C  i  0,  where  Z,  is  the  time  delay  and  p(C,)  is  the  temporal 

2  2 

normalized  covariance  of  the  field  process.  Also  pp  =  /Jpg* 

Since  (A. 5-1)  involves  an  exponential  of  an  exponential  of  an 
exponential,  and  because  a  wide  range  of  parameter  values  are  of 
interest,  care  must  be  taken  in  numerical  evaluation  of  this 
covariance  and  its  transform. 

Observe  first  that 

ko(0)  »  1  since  p(0)  *  1  .  (A. 5-2) 

Also,  as  delay  C  -♦  +<«>,  then  p  -»  0,  giving 

k^C®)  =  exp[-r^  -  2A  +  2A  exp[-p2/A)]  ^  0  .  (A. 5-3) 

The  spectrum  of  interest  is  given  by 
00 

Wq(u)  »  4  J  dC  cos(wC)  J«q(C)  for  w  i  0  ;  w  =  2nf  .  (A. 5-4) 

0 

The  nonzero  value  of  (A. 5-3)  at  C  =  “  leads  to  an  impulse  in 
spectrum  W^Cw)  at  w  «  0.  This  limiting  value,  k^(®),  must  be 
subtracted  from  covariance  (A. 5-1)  prior  to  the  numerical 
Fourier  transform  indicated  by  (A. 5-4). 
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2 

For  Mp  <<  If  the  term 


exp(-  Pp  [1  -  exp{-C)]) 

(A. 5-5) 

approaches  its  limiting  value  at  C  ®  +“  as 

follows : 

exp(-  Pp  [1  -  exp(-2.)]]  -  exp(- 

=  exp[-  p^j  [exp[r^  Pp  exp(-C) 

]-i]. 

~  exp(-  Pp)  Pp  exp(-C.) 

• 

(A. 5-6) 

This  is  a  fairly  rapid  decay  with  C  and  will  not  lead  to 
numerical  difficulty  when  ri  <<  1. 

o  Ac 

For  large  bpp/A,  the  term 

exp(-  -f  [1  -  exp(-b2;)l)  (A. 5-7) 


is  very  sharp  near  C,  ®  0;  in  factf  it  is  given  approximately  by 

2 


exp 


for  C  near  0  . 


(A. 5-8) 


Therefore,  we  define  the  sharp  component  of  covariance  ^^(2.)  as 

^2 

kg{C)  *  exp[-  A  +  A  exp|-  ~  ]  ~  exp{-A)  for  all  C,  .  (A. 5-9) 

Then 

kg(0)  »  1  -  exp(-A)  ,  kg(»)  -  0  .  (A. 5-10) 
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Now  we  let 


s  kj(C)  +  kg(C)  ,  (A. 5-11) 

where  k^(^)  is  a  flat  function  near  C,  *  0.  Then  we  can  express 
the  desired  difference  as 

ko(C)  -  k^(«»)  »  [k^{C)  -  k^C)]  +  kg(C)  - 

s  kj{C)  +  kg(C)  ,  (A. 5-12) 

where  functions  kj(C)  and  kg(C)  both  decay  to  0  at  2,  *  «.  We  now 
employ  two  separate  FFTs  on  each  of  the  functions  in  (A. 5-12). 

The  sharp  component,  k  (^),  must  be  sampled  with  a  very  small 

S 

2 

increment,  AC#  when  b/ip/A  is  large.  On  the  other  hand,  the  flat 
component 

k^CC)  -  kj(C)  -  kjjC)  (A. 5-13) 

2 

can  be  sampled  in  a  coarser  fashion.  Finally,  if  b/ip/A  is 
moderate,  we  work  directly  with  k^(J,)  -  without  breaking 

it  into  any  components. 

Two  programs  are  furnished  in  this  appendix,  one  for  moderate 
b^p/A,  and  the  other  for  the  flat  component  (A. 5-13)  when  b^p/A 
is  large.  For  sake  of  brevity,  the  Fourier  transform  of  the 
sharp  component  (A. 5-9)  is  straightforward  and  is  not  presented. 
The  particular  covariance  p(C)  adopted  is  triangular, 

p(C)  ■  1  -  for  Id  <  C,,  »  0  otherwise  ,  (A. 5-14) 

^c 

but  can  easily  be  replaced.  The  parameter  t,  is  the  cutoff  value 

w 

of  covariance  p(C)« 
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The  number  of  samples,  N,  taken  of  the  covariance,  in  order 
to  perform  the  FFT  of  (A.5-4),  is  rather  large,  so  as  to 
guarantee  a  very  small  value  of  truncation  error  at  the  upper  end 
of  the  integral,  despite  the  small  increment  AC,-  In  order  to 
keep  the  FFT  size,  Mf,  at  reasonable  values,  the  data  sequence  is 
collapsed  without  any  loss  of  accuracy  according  to  the  method 
given  in  [8;  pages  7-8]  and  [9;  pages  13  -  16].  The 
trapezoidal  rule  is  used  to  approximate  the  integral  in  (A.5-4), 
for  reasons  given  in  [8;  appendix  A]. 


10 

!  SPECTRUM  FOR  PHASE 

MODULATION  -  MODERATE 

20 

Mup=l  . 

j 

MU£ubP 

30 

Gp=. 001 

f 

Gamma'' 

40 

B£=1  . 

1 

b 

50 

A=.2 

! 

A 

60 

2c=2. *PI 

! 

Rho<2>  =  0 

for 

|2l>2c 

70 

De 1 2= . 005 

1 

Zeta  increment 

80 

M=60000 

1 

Maximum  number 

of  sarn 

90 

Mf= 16384 

! 

Size  of  FFT 

100 

DOUBLE  H.Mf.Ms.Hs 

! 

INTEGERS 

110 

DIM  X< 16384) , Y< 163 

84>  ,Cos 

<4096> 

120 

REDIM  X < 0 ; Mf - 1 ) , Y < 0 ! Mf - 1  ) 

,Cos<0;Mf/'4) 

130 

MAT  X=<0. ) 

140 

MAT  Y=<0. ) 

150 

T=2. *PI/Mf 

160 

FOR  Ms=0  TO  Mf/4 

170 

Co£<M£)=C0S<T*Ms) 

1 

QUARTER-COS 

INE 

TABLE 

180 

NEXT  Ms 

190 

T  a=Gp*Mup*Mup 

200 

IF  A=0.  THEN  220 

210 

Tb=Mup*Mup/A 

220 

Tc=2. *A*FHExp(Tb) 

230 

Ki  rif  =  FHExp<Ta+2.  -itA 

-Tc  >  ! 

CORRELATION 

AT 

INF  INI 

240 

COM  A, Bs, Zc , Ta, Tb, 

Tc , K 1 nf 

250 

T=1 . -Ki nf 

260 

PRINT  0,T 

270 

X<0)=T». 5 

1 

TRAPEZOIDAL 

RULE 

280 

FOR  Ns=l  TO  N 

290 

Corr=FNKo<N£*D6l z) 

1 

CORRELATION 

k  o  <  zet  a> 

300 

IF  N£<6  then  print  Ns,Corr 

310 

IF  ABSCCorrX  1 .  E-30  THEN 

350 

320 

M£=H£  MODULO  Mf 

1 

COLLAPSING 

330 

X<M£)=X<M£>+Corr 

340 

NEXT  Ns 

350 

PRINT  "Final  value 

of  Cor 

r  =";Corr;" 

Ns 

=  ";M£ 

360 

HAT  X=X*<D6lz*4. ) 

370 

CALL  Fft 14(Mf ,Co£<*),X<*) 

,Y<*)) 

2=26 t a 


O  f  k  O  (  2  6 1  a  ) 
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380  GINIT 

390  PLOTTER  IS  "GRAPHICS" 

400  GRAPHICS  ON 

410  UINDOW  -2, 2, -60,0 

420  LINE  TYPE  3 

430  GRID  1,10 

440  LINE  TYPE  1 

450  Delf=l./<Mf*Delz) 

460  FOR  Ms=l  TO  Mf/2 

470  F=Ms*De1f  !  FREQUENCY 

480  PLOT  LGT<F), 10.*LGTCX<Hs)) 

490  NEXT  Ms 

500  PENUP 

510  PAUSE 

520  END 

530  ! 

540  DEF  FNExp<Xtni  nus)  !  EXP<-X)  WITHOUT  UNDERFLOW 

550  IF  Xminus>708.3  THEN  RETURN  0. 

560  RETURN  EXP ( -Xm i nus ) 

570  FNEND 

580  ! 

590  DEF  FNKo<2eta)  !  CORRELATION  ko(zeta> 

600  CON  A,B£,Zc,Ta,Tb,Tc,Kinf 

610  Rho=MAX<0. , 1 . -2eta/2c )  !  TRIANGULAR  RHO 

620  El=Ta*<l.-FNExp<2eta)) 

630  E2=Tb*<l.-FNExp<Bs*2eta)) 

640  E3=A*Rho*FNExp<E2) 

650  RETURN  FNExp<El+A«<2. -Rho)-Tc»< 1 . -Rho>-E3>-Ki nf 

660  FNEND 

670  ! 

680  SUB  Fftl4<D0UBLE  N,REAL  Cos  <  ♦  >  ,  X<  » ) ,  Y<  * ) )  !  N<  =2^^  1 4=  1 6384 ;  0  SUBS 
690  !  SEE  APPENDIX  A. 3 


10 

!  SPECTRUM  FOR  PHASE  MODULATION  -  FLAT  COMPONENT 

20 

Mup=  1 . 

1  NUsubP 

30 

Gp=.O01 

!  Gamma' 

40 

Bs=l . 

1  b 

50 

A=0. 

1  R 

60 

2c=2. *PI 

!  Rho<Z)  =  0  for 

|2|>Zc5 

70 

De1z=.005 

!  2eta  increment 

80 

N=60O00 

!  Maximum  number 

of  samp 

90 

Mf=16384 

!  Size  of  FFT 

100 

DOUBLE  N,Mf,Ms,Ns 

!  INTEGERS 

110 

DIM  X<16384),Y<16384), 

Cos<4096) 

120 

REDIM  X<0:Mf-l),Y(:0iMf 

-1  ),Cos(0:Mf/'4) 

130 

MAT  X=<0.  ) 

1 46 

HAT  Y^CO./ 

150 

T»2.*PI/Mf 

160 

FOR  Ms  =  0  TO  Mf/'4 

170 

Cos<Ms)=COS<T»Ms) 

!  QUARTER-COSINE 

TABLE 

180 

NEXT  Ms 

190 

T  a*Gp*Mup«Mup 

200 

IF  A=0.  THEN  220 
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210  Tb=Mup*Mup^fl 

220  Tc=2.  *l=l*FMExp<Tb) 

230  Tb=5.E55 

240  Ki  nt'  =  FNExp(Ta+2.  *fl-Tc  )  !  CORRELftTIOM  AT  INFINITY 

250  Ea=FNExp(fi> 

260  Tbb=Tb*Bs 

270  COM  fl,Bs,2c,Ta,Tb,Tc,Kirif,Ea,Tbb 

280  T=1 . -Ki nf-< 1 . -Ea>  !  SUBTRACT  SHARP  COMPONENT 

290  PRINT  0,T 

300  X(0)=T*.5  !  TRAPEZOIDAL  RULE 

310  FOR  Ns=l  TO  N 

320  Corr=FNKl<Ns*De1z)  !  CORRELATION  LKzeta) 

330  IF  Ns<6  THEN  PRINT  Ns.Corr 

340  IF  ABSCCorrX  1 .  E-30  THEN  380 

350  Ns=N£  modulo  Mf  !  COLLAPSING 

360  X<Ms)=X<Ms)-*-Corr 

370  NEXT  Ns 

380  PRINT  "Final  value  of  Corn  =";Corr;"  Ns  =";Ns 
390  MAT  X=X*<Del2*4.  ) 

400  CALL  Fftl4(Nf,Co£<*),X<*),Y<*)) 

410  GINIT 

420  PLOTTER  IS  "GRAPHICS" 

430  GRAPHICS  ON 

440  WINDOW  -2, 2, -60,0 

450  LINE  TYPE  3 

460  GRID  1,10 

470  LINE  TYPE  1 

480  Delf=l./'<Mf*Delz) 

490  FOR  Ms=l  TO  Mf/2 

500  F=Ms*Delf  !  FREQUENCY 

510  T=X<Ms) 

520  IF  T>0.  THEN  550 

530  PENUP 

540  GOTO  560 

550  PLOT  LGT<F), 10. *LGT<T) 

560  NEXT  Ms 

570  PENUP 

580  PAUSE 

590  END 

600  ! 

610  DEF  FHExp<Xmi nus)  !  EXP(-X)  WITHOUT  UNDERFLOW 

620  IF  Xminu£>708.3  THEN  RETURN  0. 

630  RETURN  EXP<-Xminus) 

640  FNEND 

650  I 

660  DEF  FNKKZeta)  I  CORRELATION  kKzela) 

670  COM  A, Bs , Zc , Ta, Tb, Tc , K i nf , Ea, Tbb 

680  Rho=MAX<0. , 1 . -Zeta/Zc )  !  TRIANGULAR  RHO 

690  El=Ta*a.-FNExp<Zeta)> 

700  E2=Tb*< 1 . -FNExp<Bs*Zeta)) 

710  E3=A*Rho*FNExp<E2) 

720  E4=FNExp(Tbb*Zeta) 

730  Sharp=FNExp(A*< 1 . -E4) >-Ea  !  ks<zeta> 

740  RETURN  FNExp < E 1 +A* < 2 . -Rho ) -Tc * < 1 . -Rho ) -E3 > -K i nf-Shaf p 
750  FNEND 

760  ! 

770  SUE  Ffll4<D0UBLE  N,RERL  Cos < # > , X < * > , Y ( * ) >  !  N< »2" 1 4= 1 6384 ; 
780  !  SEE  APPENDIX  A. 3 


0  SUBS 
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APPENDIX  A. 6  -  EVALUATION  OF  FREQUENCY  MODULATION 
INTENSITY  SPECTRUM 

The  normalized  covariance  function  for  frequency  modulation 
is  given  by  (2.48)  in  the  main  text/  namely 

»  exp[-  ^J^  [exp(-C)  +  C 
2 

(  ^F 

+  A  p(0  exp - ^  [exp(-bC)  +  bC 

^  Ab^ 

where  C  is  the  time  delay  and  p(C)  is  the  temporal  normalized 

2  2 

covariance  of  the  field  process.  Also,  Pp  =  Ppg  and 
b  =  .  Since  (A.  6-1)  involves  an  exponential  of  an 

exponential  of  an  exponential,  and  because  a  wide  range  of 
parameter  values  are  of  interest,  care  must  be  taken  in  numerical 
evaluation  of  this  covariance  and  its  transform. 

Observe  that 

^^(0)  *  1  ,  because  p(0)  =  1  .  (A. 6-2) 

Also,  as  delay  C  then  p  0,  giving 

k^(C)  ~  exp(-  pp  (C.  -  1)  -  2a)  h  k^(C)  for  C  >  0  .  (A. 6-3) 

2 

This  term,  kj^(C)/  decays  slowly  with  Z,  if  <<  1  . 

The  spectrum  of  interest  is  given  by 


-  1]  -  A[2  -  p(C)]  + 

-  1]]]  for  2  0  ,  (A. 6-1) 


Wq(w)  *  4  J  dC  cos(wC)  for  w  i  0  ;  w  =  2nf  .  (A. 6-4) 

0 

The  spectrum  corresponding  to  the  limiting  component,  kj^(C)  in 
(A. 6-3),  is  directly  available  in  closed  form  as 


69 


TR  8887 


W^(w)  -  4  J  dC  cos(wC)  kj(C)  - 
0 

■  4 

2 

If  Kj  <<  1,  this  latter  quantity  is  large  and  very  sharply 

peaked  at  w  ■  0;  hence,  this  term  has  been  subtracted  out  and 

2 

handled  separately  when  Pp  <<  1*  The  residual  covariance, 
kQ{C)  -  k2^(Z^),  then  decays  very  rapidly  with  C,  and  is  easily 
handled  directly  by  means  of  an  FFT.  This  breakdown  is  not 
necessary  when  Pp  ~  1  and  is  avoided,  then,  by  handling 
directly  in  one  FFT. 

For  Pp/A  >>  1,  the  term 

^2 

expf - ^  [exp(-bC)  +  bC  -  I]]  (A. 6-6) 

'  Ab*^  * 

inside  the  exponential  in  (A. 6-1}  behaves  like 

^2 

exp(-  i  near  C  -  0  ,  (A. 6-7) 

where  its  major  sharp  contribution  arises.  For  example,  if 
Pp  ■  50,  A  «  1,  then  increment  AC  ■  .005  leads  to  values  for 
(A. 6-7)  of 

exp(-0.156  n^)  at  C  ■  n  AC  t  (A. 6-8) 

which  is  adequately  sampled  in  order  to  track  its  dominant 

contribution;  the  actual  sequence  of  values  is  1,  .856,  .536, 

2 

.247,  .083.  For  smaller  Pp/A,  this  sampling  interval  is  more 
than  adequate. 


>  4 


(A. 6-5) 
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Two  programs  are  furnished  in  this  appendix,  one  each  for  the 

2 

cases  of  large  and  small  The  particular  covariance  p{C) 

adopted  is  triangular, 

p(C)  =  1  -  for  Id  <  t  0  otherwise  ,  (A. 6-9) 

‘’c  ^ 

but  can  easily  be  replaced.  The  parameter  is  the  cutoff 
covariance  value  and  is  specified  numerically  in  the  examples  in 
figures  3.9a  through  3.10b. 


le 

!  SPECTRUM  FOR  FREQUENCY 

modulation  -  LA 

20 

Muf=50. 

!  MUsubF 

30 

Gp=. 001 

>  Gam  in  a' 

40 

Bs=  1 . 

!  b 

50 

ft=.  2 

!  A 

60 

2c=2. »PI 

•  Rho(Z>  =  0 

70 

Lei z=. 005 

!  Zeta  increm 

80 

M=60000 

!  Maximum  num 

90 

Mf= 16384 

!  Size  of  FFT 

100 

DOUBLE  N,Mf,Ms,Ns 

INTEGERS 

110 

DIM  X<16384),Y<16384), 

Cos<4096) 

120 

REDIM  X<0:Mf-l>,Y<0:Mf 

-l),Co£<0sMf/4> 

130 

HAT  X=<0. ) 

140 

MAT  Y=<0. ) 

150 

T  =  2.*PI/'Mf 

160 

FOR  Ms=0  TO  Mf/4 

170 

Cos<M£)=C0S<T»Ns) 

!  QUARTER-COS 

180 

NEXT  Ms 

190 

T  a=Gp*Muf  *Muf 

200 

IF  R=0.  THEN  220 

210 

Tb=Huf*Muf/'<fl*Bs*Bs) 

220 

Tc=FHExp<2. *fl-Ta)*Ta 

230 

Td=Ta*Ta 

240 

con  fl,Bs,Zc,Ta,Tb 

250 

X<0>  =  .5 

1  TRAPEZOIDAL 

260 

FOR  Ns=l  TO  N 

270 

Corr=FNKo<Ns*Del2) 

!  CORRELATION 

280 

IF  Corr<l.E-20  THEN  320 

290 

M£=N£  modulo  Mf 

!  COLLAPSING 

300 

X<Ms)=X<Ms)+Corr 

310 

NEXT  Ns 

320 

PRINT  "Final  oalue  of 

Corn  =";Corr;" 

330 

MAT  X=X*<Delz) 

340 

CALL  Ff  t  14<Mf ,  Cosi:*)  ,  X<«  )  ,  Y<*)  ) 

LRRGE  Gp^MMf''2 


for  I  2  I >2c ;  Z  =  zet  a 
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350  GIHIT 

360  PLOTTER  IS  "GRftPHICS" 

370  GRfiPHICS  OH 

380  WINDON  -4,2,-70,30 

390  LINE  TYPE  3 

400  GRID  1,10 

410  LINE  TYPE  1 

420  Del  f=l . /'<Mf*De  1  2) 

430  FOR  M£=1  to  Mf/2 

440  F  =  M£»Iielf  !  FREQUENCY 

450  T=XCM£) 

460  IF  T>0.  THEN  490 

470  PENUP 

480  GOTO  500 

490  PLOT  LGT<F) , 10. *LGT<T> 

500  NEXT  Ms 

510  PENUP 

520  fldd  =  X<0)-Tc/Td  !  ORlGi:;  CORRECTION 

530  F=l.E-4 

540  FOR  Ns=l  TO  100 

550  W=2.*PI*F 

560  Wl=Tc/'<Td  +  W*W) 

570  T=Wl+fldd 

580  IF  T>0.  THEN  610 

590  PENUP 

600  GOTO  620 

610  PLOT  LGTCF) , 10. »LGT<Wl+fldd> 

620  F*F*1.1  !  FREQUENCY 

630  IF  F>Delf  THEN  650 

640  NEXT  Ns 

650  PENUP 

660  PAUSE 

670  END 

680  ! 

690  DEF  FNExp<Xmi nus)  !  EXP<-X>  WITHOUT  UNDERFLOW 

700  IF  Xriiinus>708. 3  THEN  RETURN  0. 

710  RETURN  EXP<-Xminus) 

720  FNEND 

730  ! 

740  DEF  FNKo(2eta)  !  CORRELATION  koCzeta) 

750  COM  A,Bs,Zc,Ta,Tb 

760  E 1 =FNExp < Zei a) +Zet a- 1 . 

770  T=Bs«2eta 

780  E2=FNExp<T)+T-l . 

790  Rho=MRX<0. , 1-ZeiaxZc )  !  TRIANGULAR  RHO 

800  T=Ta*El+A*<2. -Rho)-A«Rho*FNExp<Tb*E2) 

810  RETURN  FNExp<T) 

820  FNEND 

830  ! 

840  SUB  Fftl4<D0UBLE  N,RERL  Cos < * ) , X< * ) , Y < # ) )  !  N< =2 ' 1 4  =  1 6384 5  0  SUE 
850  !  SEE  APPENDIX  A. 3 
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10 

26 

30 

40 

50 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

366 


MODULATION 
MUsubF 
Gamma" 
b 
A 

RhoCZ) 


2eta  increment 
Max  i  mum  nutuber 
Size  of  FFT 
INTEGERS 


;>,y<:&i92),Co£<:2048> 


■1  )  ,  Y<0:  Mf-1  >  , 


Cos(0: Mf /4  > 


SPECTRUM  FOR  FREOUENCY 
Muf=l  . 

Gp= . 00 1 
B£=1  . 

R=.2 

2c=2. *PI 
Del z=. 005 
N= 10000 
Mf=8192 

DOUBLE  N,Mf,M--, 

DIM  X(819 
REDIM  X<0;Mf 
MAT  X=<0.  ) 

MAT  Y=<0. ) 

T  =  2.  *PI/"Mf 
FOR  Ms=0  TO  Mf/4 
Cos<Ms)=COS<T*Ms) 

NEXT  Ms 
T  a=Gp*Muf *Muf 
IF  A=0,  THEN  220 
Tb=Muf*Mufx<A*Bs*Bs) 

T=FNExp<2. *A-Ta) 

Tc=T*Ta 
Td=Ta*Ta 

Del f=. l*Tax<2. *PI ) 

COM  A, Bs, 2c , Ta, Tb 
X<0)=. 5*< 1 . -T) 

FOR  Ns=l  TO  N 
Corr  =  FNKol':Ns*Delz) 

IF  ABSCCorrX  1 .  E-30  THEN 
Ms=Ns  MODULO  Mf 
X(Ms)=X<Ms)+Corr 
NEXT  Ns 

PRINT  "Final  value  of  Corn  =";Corr; 
NAT  X=X*<Delz) 

CALL  Ffl 14<Nf ,Co£<*),X<*),Y<#)> 


SMALL  Gp*Nuf''2 


=  0  for 


|2l>2c;  Z=zeta 
of  sample: 


f  kc 


!  QUARTER-COSINE  TABLE 


!  INCREMENT  IN  FREQUENCY 


!  TRAPEZOIDAL  RULE 


!  CORRELATION 
340 

!  COLLAPSING 


ko<zeta)-kl<zeta> 


Ns 


Ns 


I.'  z  e  t  a  > 
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370  GINIT 

380  PLOTTER  IS  "GRAPHICS" 

390  GRAPHICS  OH 

400  WINDOW  -4,2,-70,30 

410  LINE  TYPE  3 

420  GRID  1,10 

430  LINE  TYPE  1 

440  FOR  M£=1  to  2000 

450  F  =  M£.*Delf  !  FREQUENCY 

460  W=2.*PI*F 

470  Wl=Tc/<Td+W*W)  !  SHARP  SPECTRAL  COMPONENT 

480  T=Mf*De1z*F 

490  N£=INT<T) 

500  Fr=T-N£ 

510  W2  =  Fr-*X<N£+n  +  <  1  . -Fr)*X<N£)  !  BROAD  SPECTRAL  COMPONENT 

520  PLOT  LGTCF) , 10. *LGT<W1+W2  ■ 

530  NEXT  Ms 

540  Ns=MAX<N£,1) 

550  FOR  M3=N£  to  Mfv2 

560  F  =  Ms/’<Hf*D€  1 2)  !  FREQUENCY 

570  W  =  2.*PHrF 

580  Wl=Tc/'<Td  +  W*W) 

590  W2=X<N£) 

600  T=W1+W2 

610  IF  T>0.  THEN  640 

620  PENUP 

630  GOTO  650 

640  PLOT  LGT<F) , 10. *LGT<T) 

650  NEXT  Ns 

660  PENUP 

670  PAUSE 

680  END 

690  ! 

700  DEF  FNExpCXmi nus)  !  EXP(-X>  WITHOUT  UNDERFLOW 

710  IF  Xtr.irius>708.3  THEN  RETURN  0. 

720  RETURN  EXP(-Xminus) 

730  FNEND 

740  - 

750  DEF  FNKoUZeta)  !  CORRELATION  k o ( zct a) -k  1  < zet a> 

760  COM  A, Bs, 2c , Ta, Tb 

770  E 1 =FNExp < 2el a) +Zet a- 1 . 

780  T=Bs*2eta 

790  E2  =  FNExp<T:-  +  T-1  . 

800  kho=MAX<0. , l-2eta/Zc )  !  TRIANGULAR  RHO 

810  T=Ta*El+A*<2. -Rho)-A*Rho*FNExp<Tb*E2) 

820  RETURN  FNExp  <  T  > -FNExp  <  Ta*  <  Zet  a- 1 .  ) +  2 .  > 

FNEND 
840  ! 

850  SUB  Ff  i  MCDOUBLE  N,REAL  Cos  <  *  j  ,  X  <>  >  ,  Y  C  *  )  )  !  N<  =2"  1  4=  1  6384 

860  !  SEE  APPENDIX  A. 3 
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